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on  spray  have  been  studied  and  graphically  presented.  The  calculated  spray 
penetration  with  and  without  crossflow  of  air  has  been  compared  with  the 
available  experimental  data  of  other  researchers  and  good  agreement  between 
the  two  is  noticed. 

The  model  can  predict  the  rate  of  combustible  mixture  formation,  the 
rate  of  heat  release  and  cylinder  pressure  as  a  function  of  time  in  direct 
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unusable,  so  comparisons  between  modeled  and  measured  cylinder  pressure 
rise  data  are  not  included  in  this  report.  As  soon  as  such  data  is 
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NOMENCLATURE 


Surface  area  of  droplet 
Coefficient  of  discharge 
Coefficient  of  drag  for  fuel  droplet 
Nozzle  diameter 

Sauter  mean  diameter  of  droplet 

Mass  averaged  diameter  of  droplet 

Gravitational  constant 

Evaporation  rate  constant  in  stagnant  air 

Evaporation  rate  constant 

Length  of  nozzle  hole 

Mass  of  liquid  fuel 

Number  of  droplets 

Chamber  gas  pressure 

Cylinder  charge  pressure 

Injection  pressure 

Discharge  from  nozzle 

Reynolds  number 

Schmidt  number 

Time  from  start  of  injection 

Chamber  gas  temperature 

Injection  timing 

Liquid  fuel  velocity 

Jet  velocity  at  nozzle  exit 


Units 


cm/sec 
cm2  /sec 

O 

cm  /sec 


mm  /  stroke 


deg  BTDC 

cm/sec 

cm/sec 


*  Indicates  units  used  in  the  text  unless  stated  otherwise 


I 


t 


IV 

|i!' 


r 


:  Relative  velocity  of  fuel  droplet  (wrt  air  velocity)  cm/sec 


or  ‘ 

Relative  velocity  of  droplet  at  the  beginning  of 
any  interval  of  time 

cm/sec 

Penetration 

cm 

P  : 

Pressure  differential  across  nozzle 

atm 

^  : 

Chamber  gas  density 

gm/cc 

P1  : 

Liquid  fuel  density 

gm/cc 

e  : 

Spray  cone  angle 

degrees 

Tl : 

Ignition  delay  period 

seconds 

*  : 

Viscosity  of  chamber  gas 

gm/cm-sec 

4i  : 

Viscosity  of  liquid  fuel 

gm/cm-sec 

’a  : 

Kinematic  viscosity  of  chamber  gas 

2  , 

cm  /sec 

DEFINITIONS 

ds  = 

£ni  d, 

d,  -  (£"|  d,3  )°'33 

Re  =  Reynolds  number 

=  velocity  of  droplet  x  droplet  diameter 
kinematic  viscosity  of  ambient  gas 

Sc  =  Schmidt  number 

=  kinematic  viscosity  of  ambient  gas 
diffusivity  coefficient 


1 .  INTRODUCTION 


Diesel  Engines  are  used  for  a  wide  variety  of  applications,  not  only 
for  the  purpose  of  transport  but  also  as  a  source  of  power  for  stationary 
systems.  Their  wide  acceptance  as  prime  movers  is  attributed  to  their 
superior  fuel  consumption  characteristics  in  heavy  duty  vehicles,  and  low 
hydrocarbon  and  carbon  monoxide  emissions  in  lighter  duty  vehicle 
applications  although  they  exhibit  somewhat  higher  nitric  oxide  and 
particulate  emission.  The  purpose  of  this  research  project  was  to  create  a 
computer  model  of  the  spray  combustion  process  in  a  direct  injection  diesel 
engine  which  could  then  be  used  to  study  the  parameters  necessary  to 
optimize  engine  performance. 

Computer  modeling  of  the  combustion  process  is  a  valuable  complement 
to  engine  design  and  testing  since  it  allows  the  investigation  of  varying 
engine  design  parameters  without  the  expense  of  building  and  testing 
numerous  possible  design  alternatives.  The  model  developed  here  takes  a 
phenomenological  approach  to  combustion  analysis  and  relies  on  previous 
results  of  experiments  on  various  fundamental  aspects  of  the  combustion 
process.  A  fully  detailed  simulation  of  engine  operation  is  not  yet 
possible  due  to  a  lack  of  understanding  of  all  the  complex  phenomena 
occurring  in  an  operating  engine. 

To  achieve  improvements  in  the  performance  of  Direct  Injection  (D.I.) 
Diesel  Engines,  the  process  of  combustion  needs  to  be  better  understood. 

The  complicated  mechanisms  involved  in  diesel  combustion  are  the  subject  of 
on-going  research  by  many  individuals.  From  what  is  known,  there  are 
strong  indications  that  the  diesel  spray  and  mechanics  of  its  formation 
preceeding  combustion  significantly  influence  the  entire  process  and  thus 


calls  for  primary  attention  in  any  investigation. 

Fuel  spray  injection  from  a  nozzle  into  a  combustion  chamber  is 
followed  by  the  division  of  liquid  fuel  into  a  large  number  of  fine 
droplets  to  ensure  good  distribution  of  fuel  in  the  chamber  and  to  increase 
the  surface  area  through  which  rapid  heat  transfer  rates  may  be 
accomplished.  Atomization,  phenomenon  responsible  for  fuel  breakup,  is 
a  direct  result  of  any  one  or  combination  of  the  following:  cavitation, 
turbulence,  roughness  of  nozzle,  pressure  variation  across  the  nozzle  and 
the  formation  of  surface  waves.  The  emerging  droplets  come  into 
contact  with  hot  swirling  air  above  the  piston.  They  exchange  and  conserve 
mass,  momentum  and  energy  while  traveling  away  from  the  nozzle  tip, 
overcoming  drag  and  vaporizing  to  form  a  well  distrubuted  fuel  air  mixture 
in  the  combustion  chamber.  Combustion  then  takes  place  in  three 
recognizable  stages.  The  first  stage,  at  the  beginning  of  autoignition, 
corresponds  to  the  process  of  burning  of  premixed  fuel  and  air  existing 
within  specific  flammability  limits  at  the  end  of  ignition  delay.  A  rapid 
rate  of  pressure  rise  is  the  result  of  this  process.  The  second  stage  is 
the  burning  of  the  combustible  mixture  formed  during  the  rapid  pressure 
rise  period.  Finally,  in  the  third  stage,  fuel  still  in  the  liquid  phase 
and  a  very  rich  mixture  of  fuel  and  air  at  the  core  of  the  spray,  burns  in 
a  diffusion  flame. 

In  the  present  study,  efforts  were  focused  on  developing  a  computer 
model  to  predict  spray  character istics  during  its  formation,  locate  the 
regions  where  autoignition  commences  and  then  analyze  combustion  in  the 
stages  described  above  for  environmental  conditions  obtaining  in  actual 
engines.  The  computer  model  is  capable  of  predicting  the  rate  of  heat 
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release  and  cylinder  pressure  with  respect  to  the  engine  crank  shaft 
rotation.  The  model  predictions  with  and  without  crossflow  of  air 
for  spray  penetration  were  verified  from  available  experimental  data  of 
other  researchers.  Work  is  presently  under  way  to  conduct  experiments 
on  a  Ricardo  -  Hydra  Direct  Injection  Diesel  Research  Engine  procured 
by  the  department  to  validate  the  results  of  the  model. 


2.  LITERATURE  REVIEW 

An  extensive  survey  of  the  literature  was  done  to  review  previous 
work  on  studying  the  various  processes  occuring  within  the  cylinder 
and  other  modeling  techniques  adopted.  Beginning  with  Schwietzer  in  1938, 
efforts  have  been  directed  at  understanding  the  complexities  of  inter¬ 
actions  between  processes  that  follow  the  injection  of  fuel  in  diesel 
engines.  To  form  a  reasonably  strong  basis  for  any  further  investigation, 
it  was  necessary  to  review  the  opinions  of  other  researchers  regarding 
the  nature  of  spray  penetration,  effect  of  various  engine  parameters  on 
penetration,  mixing  of  fuel  vapor  in  air  and  finally  the  process  of 
combustion . 

2.1  NATURE  AND  GEOMETRY  OF  SPRAY 

Work  done  to  date  can  broadly  be  classified  as  either  analytical  or 
experimental  or  a  combination  of  both.  The  analytical  approach  adopted  by 
Ogasawara  and  Sami[1]«  and  Elkotb  and  Rafat  [2]  treated  the  spray  as  an 
agglomerate  of  droplets  with  their  studies  confined  to  one  such  droplet. 

On  the  other  hand,  Oz  [3]  and  Melton  [4]  opted  in  favor  of  a  "Mixing  Jet 
Theory"  by  regarding  the  spray  as  fuel  vapor  injected  at  high  pressure  into 
lower  pressure  air.  Experiments  conducted  on  diesel  engines  have 
revealed  the  coexistence  of  liquid  and  vapor  phases  in  actual  fuel  sprays. 
In  fact  Hiroyasu  and  Kadota  [5]  successfully  demonstrated  that  fuel 
droplets  never  attain  critical  temperature  even  when  injected  into  a 
supercritical  environment,  thereby  confirming  the  heterogeneous  character 
of  fuel  sprays  inside  diesel  engines. 

•Numbers  in  brackets  denote  the  reference  found  at  the  end  of  this 
repor t . 


There  appears  to  be  no  general  consensus  among  researchers  regarding 
the  geometry  of  spray,  especially  in  the  presence  of  swirling  air.  Chiu, 
Shahed,  and  Lyn  [6]  of  Cummins  Engine  Company  have  assumed  the  shape  to  be 
like  that  of  a  blunt  body  whereas  Hiroyasu  and  Kadota  [7]  have  considered 
the  spray  as  a  right  circular  cone  under  no  swirl  conditions. 

2.2  AIR  MOTION 

Air  motion  in  a  high  speed  direct  injection  diesel  engine  cylinder 
consists  of  swirl  motion,  an  inward  radial  component  of  air  motion, 
commonly  called  "Squish",  and  an  axial  component  of  air  velocity  due  to  the 
motion  of  the  piston.  The  role  of  moving  air  is  to  facilitate  proper 
mixing  of  the  fuel  vapor  with  air  and  enhance  the  subsequent  process  of 
combustion . 

The  inlet  valve  geometry  and  the  extent  of  its  opening  affect  the 
generation  of  swirl  in  a  combustion  chamber.  For  the  purpose  of  analysis, 
swirl  can  be  assumed  as  solid  body  rotation  about  the  cylinder  axis  [8]. 
Squish  tends  to  cause  swirl  to  deviate  from  true  solid  body  rotation. 

Squish  has  a  maximum  value  close  to  TDC  and  depends  largely  on  the  shape 
of  the  combustion  chamber. 

2.3  SPRAY  DEVELOPMENT 

Schweitzer  [9]  was  among  the  first  to  investigate  the  penetration  of 
sprays.  Using  oil  in  his  experiments  conducted  in  a  pressure  chamber  with 
stagnant  air,  he  determined  the  factors  that  have  considerable  impact  on 
spray  tip  penetration.  Dimensional  analysis  served  to  establish  a  relation 
between  penetration  and  injection  pressure,  chamber  pressure,  orifice  size 
and  shape,  and  viscosity.  It  was  a  good  beginning,  considering  that 
sophisticated  photography  equipment  was  not  available  in  1938,  and  the  fact 


that  he  ignored  many  important  factors  such  as  evaporation,  effect  of 
elevated  temperatures  and  air  swirl.  Schweitzer’s  correlation,  along  with 
those  of  others  discussed  below,  have  been  summarized  in  tabular  form  in 
Table  1 . 

Since  the  1960's  a  lot  of  research  has  been  conducted  on  this  topic. 
Parks,  et  al  [10]  performed  experiments  which  helped  them  conclude  that 
higher  gas  temperature  reduces  the  penetration  of  spray  because  of  rapid 
evaporation,  provided  the  chamber  density  remains  constant.  They  also 
found  that  a  larger  orifice  diameter  increases  penetration  due  to  increased 
droplet  size.  The  effects  of  varying  injection  pressure  and  chamber 
pressure  were  also  included  in  the  scope  of  their  investigation.  While 
higher  injection  pressure  increases  the  velocity  of  droplets  which  tends 
to  increase  the  penetration,  it  also  results  in  the  formation  of  smaller 
droplets  thereby  reducing  the  time  of  evaporation  and  consequently  the 
penetration  of  droplets.  The  net  effect  is  a  small  Increase  in  the  value 
of  penetration.  Similarly,  a  lower  chamber  pressure  increases  the  pressure 
differential  across  the  nozzle,  producing  smaller  droplets  and  tending  to 
decrease  penetration.  But,  as  drag  is  also  decreased  the  final  result  is  a 
slight  increase  in  penetration.  By  curve  fitting  their  experimental  data, 
they  modified  Schweitzer's  correlation.  However,  they  failed  to  recognize 
the  importance  of  momentum  exchange  between  fuel  and  surrounding  air. 

Ogasawara  and  Sami  [l]  proposed  an  expression  to  evaluate  spray  tip 
penetration  after  studying  the  behavior  of  a  single  droplet.  While 
viscous  drag  and  exchange  of  momentum  with  the  surrounding  air  was  kept 
in  mind,  the  effect  of  any  evaporation  of  the  droplet  during  motion  in  a 
high  temperature  and  pressure  environment  was  assumed  to  be  negligible. 


Oz  [3]  in  1969  was  among  the  first  few  to  advocate  the  concept  of 
analyzing  fuel  sprays  in  diesel  engines  as  being  equivalent  to  injecting 
vapors  at  high  speed  into  air.  The  empirical  formula  that  Oz  presented 
gave  a  fairly  good  match  with  experimental  values  for  low  back  pressures. 

Melton  [4]  proposed  that  at  a  distance  from  the  source  of  about  30 
nozzle  diameters,  a  diesel  spray  is  essentially  an  air  jet  bearing  a  fog 
of  fuel  droplets.  As  such,  Melton  applied  the  "Theory  of  Submerged  Jets" 
to  spray  mechanics  and  derived  an  expression  for  spray  penetration. 

In  1971,  Dent  [11]  reviewed  all  the  experimental  work  done  in  cold  and 
hot  bombs  and  utilized  the  information  to  generate  an  equation  for  fuel 
penetration  under  cold  bomb,  hot  bomb,  and  engine  working  conditions. 

Later  Hay  and  Jones  [12]  recommended  the  use  of  this  equation  for  all 
conditions  except  for  very  high  gas  densities. 

An  attempt  to  model  fuel  spray  trajectory  under  conditions  pertinent 
to  actual  engine  cylinders  was  made  in  1977  by  Elkotb  and  Rafat  [2].  The 
effect  of  shape  of  the  combustion  chamber,  air  velocity  pattern,  heat 
transfer  and  drag  force  variations  were  all  considered  in  the  model 
to  predict  the  location  of  a  fuel  droplet  at  any  Instant.  Using  basic 
equations  of  mass,  momentum  and  energy  conservation,  non-linear  differen¬ 
tial  equations  were  developed.  These  were  solved  numerically  and  the 
results  were  verified  by  experimental  investigation.  Their  model  sought 
to  improve  upon  the  ballistic  approach  adopted  by  Ogosowara  in  1966. 

All  experimental  work  done  before  1978  gave  poor  results  for  the 
initial  port  of  the  injection.  But  it  is  this  initial  period  of  about  one 
millisecond  (MS)  that  is  of  interest  in  an  actual  diesel  engine  since  it 
corresponds  closely  to  the  ignition  delay  period.  Hiroyasu,  et  al  [13] 


developed  techniques  to  overcome  this  drawback.  They  performed  experiments 


with  quiescent  air  as  well  as  with  swirling  air.  Study  of  fuel  penetration 
was  done  under  varying  ambient  gas  pressure,  injection  pressure,  nozzle 
geometry  and  ambient  gas  temperature.  It  was  concluded  that  spray 
penetration  is  proportional  to  time  in  the  early  stage  of  injection  and 
thereafter  to  the  square  root  of  time.  Thus  the  fuel  jet  was  divided  into 
a  developing  region  and  a  fully  developed  region,  with  transition  occurring 
at  a  time  called  the  breakup  time. 

Photographs  taken  by  Kuniyoshi,  et  al  [14]  in  1980  revealed  that  the 
spray  is  composed  of  two  distinct  zones,  namely  the  main  jet  region  and  a 
mixing  region.  The  main  jet  region  lies  towards  the  core  of  the  spray  and 
holds  small  sized  droplets  moving  at  high  velocities,  while  the  mixing 
region,  containing  somewhat  larger  droplets  moving  with  lower  velocity, 
encircles  the  main  jet  region.  They  observed  that  the  total  penetration 
could  be  subdivided  into  an  initial  part,  a  mixing  part  and  a  dilution 
portion.  Further,  for  high  charge  pressure  and  varying  elevated 
temperatures  the  initial  part  of  the  penetration  remained  unchanged. 

A  new  dimension  was  added  by  Varde,  Popa  and  Varde  [15]  in  the  search 
for  independent  parameters  that  influence  penetration.  They  established 
the  dependence  of  penetration  on  nozzle  L/D  ratio  in  addition  to  the 
indirect  effect  it  has  on  altering  the  coefficient  of  discharge.  Working 
at  extremely  high  injection  pressures,  they  also  found  that  the  decay  in 
velocity  was  directly  proportional  to  the  injection  pressure  for  high 
penetration  values  but  showed  a  reversed  trend  for  smaller  values  of 
penetration . 

A  study  of  the  process  of  atomization  using  diesel  fuel  injection  at 
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high  pressure  and  high  temperature  into  a  quiescent  gaseous  environment  was 


undertaken  by  Takeuchi,  et  al  [16]  in  1383.  They  explained  that  fuel,  when 
injected  from  any  nozzle,  has  a  jet  cone  which  grows  in  a  wavy  form  then 
breaks  up  into  thin  threads.  Finally  small  droplets  are  formed  from 
these  filaments.  The  spray  tip  penetration  was  seen  by  them  to  be  composed 
of  3  linear  relations.  The  first  corresponded  to  the  period  required  for 
jet  breakup,  the  second  was  present  till  the  disappearance  of  wavy  flow  and 
the  third  gave  the  spray  tip  penetration  beyond  that  period.  All  three 
equations,  one  for  each  phase,  are  given  in  Table  1. 

A  "Stochastic  Thick  Spray  Model"  has  been  developed  in  the  form  of  a 
computer  code  by  Kuo,  Yu,  and  Shahed  [17]  to  quantify  fuel  and  air  in 
different  regions  in  space.  The  analysis  treats  the  spray  os  being 
composed  of  discrete  computational  particles  representing  a  group  of 
droplets  and  follows  them  along  lagrangian  trajectories  accounting  for 
mass,  momentum  and  energy  exchange  at  every  time  step. 

Fragoulis  and  Henein,  in  a  recent  publication  [18],  have  examined  most 
of  the  existing  expressions  for  spray  penetration  and  derived  a  correlation 
of  their  own  by  "averaging"  the  powers  to  which  independent  parameters  have 
been  raised  in  these  expressions.  They  also  introduced  a  new  concept  of 
"Mean  Penetration  Diameter"  allowing  description  of  spray  tip  penetration 
by  knowing  the  initial  diameter  of  the  droplets. 

Work  is  under  progress  at  the  University  of  Manchester  England  [19]  to 
examine  the  formation  of  spray  in  stagnant  air  as  well  as  air  moving  in  a 
direction  perpendicular  to  the  direction  of  injection.  The  latter  has  been 
achieved  in  a  wind  tunnel  with  induced  air  velocities  of  8  and  12  m/s.  The 
experimental  data  has  shown  spray  penetration  to  be  mainly  dependent  on  gas 
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density  with  a  relatively  weak  dependence  on  gas  temperature.  At  room 
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temperature,  no  substantial  difference  was  recorded  in  the  values  with  and 
without  crossflow. 

2.4  IGNITION  DELAY  AND  RATE  OF  BURNING 

Although  numerous  definitions  for  ignition  delay  have  been  given  by 
different  authors,  all  seem  to  agree  that  it  is  the  period  extending  from 
the  beginning  of  injection  to  the  time  of  measurable  combustion.  Henein 
and  Bolt  [20]  found  that  pressure  rise  was  the  first  noticeable  change 
occurring  in  diesel  combustion  and  emphasized  the  use  of  measuring  ignition 
delay  in  terms  of  pressure  rise.  Their  engine  tests  demonstrated  that 
increase  in  cylinder  air  pressure,  fuel/air  ratio,  cooling  water  tempera¬ 
ture  and  engine  speed,  all  work  towards  shortening  the  delay  period. 

The  same  authors,  in  1969  [21],  suggested  an  exponential  correlation 
for  evaluating  ignition  delay  in  terms  of  apparent  activation  energy  and 
mean  integrated  temperature  of  air  charge.  An  increase  in  activation 
energy  as  per  that  correlation,  increases  ignition  delay  but  any  increase 
in  temperature  decreases  ignition  delay. 

Spontaneous  ignition  delay  of  fuel  droplets  and  fuel  spray  in  a  high 
pressure  gaseous  environment  was  examined  closely  by  Hiroyasu  and  Kadota  in 
two  separate  experiments  [5].  The  effect  of  charge  pressure,  temperature 
and  oxygen  concentration  of  the  ambient  gas  on  ignition  delay  was  recorded 
and  used  to  yield  an  expression  for  finding  ignition  delay.  The  constants 
in  the  expression  were  found  to  be  functions  of  fuel  properties  and  their 
values  have  been  determined  for  various  fuels.  Similar  effects  of  varying 
intake  air  pressure  and  temperature,  cooling  water  temperature  and  engine 
speed  on  ignition  delay  in  a  single  cylinder,  direct  injection,  4  stroke 


diesel  engine  were  observed  by  Wong  and  Steere  in  1982  [22]. 

A  number  of  ignition  delay  models  for  heterogeneous  charge  engines 
were  compared  by  Primus  and  Wong  [25]  in  1985.  They  observed  that  there  is 
widespread  agreement  that  the  physical  delay  is  very  short  relative  to  the 
chemical  delay;  and,  therefore,  the  total  delay  may  be  measured  using  the 
Arrhenius  equation.  In  the  majority  of  cases,  mean  cylinder  pressure  and 
temperature  are  used  for  computing  the  ignition  delay  period.  As  far  as 
heat  release  is  concerned,  many  of  the  existing  models  consider  an  apparent 
heat  reler se  diagram  as  a  part  of  their  program  input  to  calculate  the 
actual  heat  release  by  iterating  to  satisfy  the  burning  rote  equation  and 
the  equations  of  conservation  of  energy  simultaneously. 
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3.  DIESEL  ENGINE  COMBUSTION  CHARACTERISTICS 


Inspite  of  the  complexities  involved  in  diesel  combustion,  a  large 
number  of  correlations  and  models  have  been  proposed  in  the  past,  but  very 
seldom  do  we  find  them  being  applied  under  actual  engine  conditions.  Most 
theoretical  work  assumes  the  properties  of  chamber  gas  and  fuel  to  be  a 
constant  and  experiments  have  invariably  been  performed  in  cold  or  hot 
chambers,  having  constant  volume,  fixed  temperature  and  non  varying 
pressure  drop  across  the  nozzle  for  each  test  run.  Furthermore,  the 
photographic  technique  employed  in  all  experiments  gives  good  results  for 
time  greater  than  1  ms  from  the  start  of  injection.  As  mentioned  earlier, 
it  is  the  time  prior  to  1  ms  that  is  of  significance  since  high  speed 
diesel  engines  rarely  display  an  ignition  delay  of  over  1  ms.  Thus,  the 
relationships  are  seldom  used  in  determining  spray  characteristics  during 
time  less  than  one  millisecond  from  the  start  of  injection. 

The  expressions  for  spray  penetration  summarized  in  Table  1  cannot  be 
compared  directly,  as  each  was  derived  under  a  different  set  of  environment 
conditions.  Nevertheless,  a  close  examination  of  each  expression  provides 
insights  about  the  effects  of  various  parameters  on  spray  penetration. 

These  parameters  are  listed  on  the  next  page  and  a  brief  discussion  on  each 
follows  in  the  next  section. 

AP  Mean  effective  pressure  across  the  nozzle 

:  Coefficient  of  discharge  for  the  nozzle 

d  :  Diameter  to  nozzle  opening 

L/d  :  Length  to  diameter  ratio  of  the  nozzle 

p  Ambient  gas  density 

cl 

T 

a 


Ambient  gas  temperature 


'M 


Liquid  fuel  density 


SMD  Sduter  mean  diameter  of  droplets 


Spray  cone  angle 


K  :  Evaporation  rate  constant 


t  :  Time  from  the  start  of  injection 

In  addition  to  the  variables  listed  above,  penetration  in  an  actual 
diesel  engine  may  also  be  affected  by  engine  dimensions,  engine  speed, 
duration  of  injection,  swirl  ratio,  angle  of  inclination  of  nozzle, 
distribution  of  SMD  in  spray  cone,  distribution  of  fuel  by  mass  across  the 
cone  and  some  fuel  properties. 


3.1  EFFECT  OF  SOME  PARAMETERS  ON  PENETRATION 


(1)  Pressure  differential  across  nozzle 


The  pressure  differential  across  a  nozzle  can  be  altered  either 
by  changing  the  injection  pressure  or  by  varying  the  chamber  pressure. 

(a)  Changing  injection  pressure:  An  increase  in  injection 
pressure  increases  the  velocity  of  droplets  which  tends  to  increase 
the  penetration.  At  the  same  time,  higher  injection  pressure  also 
results  in  the  formation  of  smaller  droplets  having  shorter  life  and 
hence  less  penetration.  The  net  effect,  as  noticed  by  many,  is  a 


slight  increase  in  the  value  of  penetration. 

(b)  Altering  chamber  pressure:  Higher  chamber  pressure  results 
in  the  formation  of  larger  droplets  capable  of  travelling  farther,  bu 


since  the  drag  on  the  droplets  also  increases,  the  net  result  is  a 
slight  decrease  in  the  penetration. 

(2)  Coefficient  of  discharge  for  the  nozzle 

Coefficient  of  discharge  has  a  very  obvious  effect  on 


»c^n,r 


penetration.  Penetration  increases  as  coefficient  of  discharge  is 
increased . 

(3)  Orifice  diameter 

With  an  increase  in  orifice  diameter,  larger  sized  droplets  are 
formed.  They  can  survive  longer,  and  even  though  the  drag  they  experience 
escalates,  they  move  much  further  as  compared  to  smaller  droplets. 

(4)  L/d  ratio  of  the  nozzle 

A  higher  L/d  ratio,  in  addition  to  increasing  the  coefficient  of 
discharge  for  the  nozzle,  also  inhibits  the  decay  in  center  line  velocity 
thereby  increasing  the  penetration. 

(5)  Ambient  gas  density 

For  all  other  parameters  remaining  unchanged,  an  increase  in  chamber 
gas  density  increases  the  drag  on  fuel  droplets,  resulting  in  smaller 


penetration . 

(6)  Ambient  gas  temperature 

Although  an  increase  in  the  ambient  gas  temperature  decreases  the  life 
of  fuel  droplets,  the  ambient  gas  density  is  also  decreased  and 
subsequently  reduced  drag  enhances  penetration. 

The  remaining  parameters  affecting  penetration  are  dependent  variables 
and  can  be  expressed  in  terms  of  the  above  discussed  parameters.  Knowing 
their  relationship  to  the  above  six,  their  effects  can  also  be  predicted. 

The  variables  mentioned  in  (1)  through  (6)  above  are  characteristics 
of  the  engine  or  the  working  environment.  Their  value  serves  as  an 
essential  input  to  any  model  that  portrays  the  development  of  fuel  sprays. 


17 


3.2  AN  INSIGHT  TO  PROBABLE  INPUTS 


Souter  Mean  Diameter  (SMD),  the  spray  cone  angle  (9  ),  the  evaporation 
rate  constant  (K),  and  the  ignition  delay  period  (t )  are  the  few  other 
possible  inputs  that  need  to  be  scrutinized.  Vital  information  is 
available  on  them  and  deserves  individual  inspection  in  the  following 
paragraphs . 

3.2.1  SAUTER  MEAN  DIAMETER 


Several  mean  diameters  of  droplets  in  a  spray  have  been  used  to 
describe  the  quality  of  atomization,  but  Sauter  Mean  Diameter  is  of 
particular  interest  in  combustion  application.  The  definition  of  SMD 
relates  the  surface  area  to  the  volume  of  the  droplets,  thereby  providing 
a  very  essential  relationship  between  surface  area  and  volume  that  is 
required  for  analyzing  evaporation. 

The  variation  of  SMD  with  location  in  a  spray,  and  the  effect  of 
injection  pressure,  chamber  pressure  and  nozzle  geometry  on  this  variation 
have  been  studied  by  Hiroyasu  and  Kadota  [26],  The  empirical  formula 


proposed  by  them,  and  commonly  accepted  for  use,  is  given  by 


SMD  =  C  (Ap) 


-0.1  35 


,0.1  21 


,0.1  31 


(3.1  ) 


where, 


SMD  =  Sauter  Mean  Diamter  in  pm 
AP  =  Mean  effective  pressure  drop  in  MPa 

p  =  Air  or  chamber  gas  density  in  Kg/m3 

cl 

Q  =  Amount  of  fuel  delivered  in  mm3 /stroke 
C  =  Constant  depending  on  nozzle  type 
*=  25.1  for  pintle  nozzle 
=  23.9  for  hole  nozzle 
»  22.4  for  throttling  nozzle 


A  typical  distribution  of  the  SMD  across  the  spray  as  obtained  from 


the  above  reference  is  shown  in  Figure  1 . 


3.2.2  SPRAY  CONE  ANGLE 


Bracco  [27]  has  shown  that  in  the  case  of  high  speed  diesel  engines 
the  atomization  of  liquid  fuel  is  spontaneous  and  the  breakup  length  is 
practically  non-existent.  In  other  words,  divergence  of  spray  begins 
immediately  at  the  exit  of  the  nozzle.  The  same  observation  has,  however, 
not  been  reported  for  low  speed  diesel  engines  or  for  injection  into 
constant  volume  hot  and  cold  chambers .  The  angle  of  divergence  is  also 
known  as  the  spray  cone  solid  angle.  This  angle  is  a  function  of  nozzle 
dimensions,  chamber  gas  density  and  the  pressure  differential  across  the 
nozzle  [15] . 

Varde,  Popa  and  Varde  [15]  recently  studied  the  effect  of  chamber  gas 
density  and  AP  on  the  spray  cone  angle.  Their  results  are  represented 
graphically  in  Figure  2.  A  significant  change  in  the  behavior  of  the 
curves  is  noticed  as  the  nozzle  L/d  ratio  is  increased.  The  transition  in 
the  trend  was  found  to  occur  at  an  L/d  value  of  about  4. 

3.2.3  EVAPORATION  RATE  CONSTANT 

The  evaporation  rate  constant  for  the  liquid  droplets  of  fuel  is  a 
function  of  the  pressure  and  temperature  of  the  medium  surrounding  the 
droplets.  As  expected,  the  effect  of  higher  temperature  is  to  increase  the 
value  of  the  evaporation  rate  constant  and  consequently  decrease  the  life 
span  of  a  droplet.  At  higher  temperatures,  an  increase  in  gas  pressure 
further  increases  the  evaporation  rate,  whereas  Just  the  opposite  happens 
if  the  gas  temperature  is  low.  The  dependence  of  the  evaporation  rate  on 
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temperature  and  pressure  for  various  types  of  fuel  in  a  quiescent  air 


environment  has  been  determined  by  Hiroyasu  and  Kadota  [5]  and  is 
reproduced  graphically  in  Figure  3. 

3.2.4  IGNITION  DELAY  PERIOD 

Depending  on  the  definition  of  ignition  delay  various  correlations  for 
computing  the  delay  period  have  been  proposed.  Since  it  is  widely  accepted 
to  consider  the  beginning  of  rapid  pressure  rise  to  be  the  end  of  ignition 
delay,  the  work  done  by  Hiroyasu  and  Kadota  [5]  is  of  particular  interest. 
They  have  conducted  experiments  with  liquid  sprays  and  single  droplets  at 
high  temperature  and  pressure  to  establish  the  degree  of  influence  of 
chamber  gas  pressure,  gas  temperature  and  oxygen  concentration  in  the 
ambient  gas  on  ignition  delay.  The  following  equation  is  the  outcome  of 
their  efforts: 

t  =  A  (P  )B  (<t)C  exp  (D/T  )  (3.2) 

e  a 

where , 

Pa  =  Chamber  gas  pressure 

=  Oxygen  concentration  in  ambient  gas 
=  Mean  temperature  of  gas 

A,  B,  C,  and  D  are  constants,  whose  values  for  some  commonly  available 
fuels  are  also  given  in  reference  5. 
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4.  MODEL  DEVELOPMENT 


4.1  MODEL  ASSUMPTIONS 

The  following  assumptions  were  made  while  developing  the  computer 
model  to  predict  the  spray  geometry,  the  fuel-air  distribution,  the  time  of 
ignition  delay  and  the  subsequent  pressure  rise  during  combustion: 

(1)  The  diesel  fuel  spray  has  been  assumed  to  be  heterogeneous  in 
character,  with  spray  disintegration  and  droplet  formation  occurring  at 
the  nozzle  exit.  This  assumption  is  consistent  with  the  findings  of 
Bracco  [27]  for  high  speed  diesel  engines. 

(2)  The  mean  value  of  the  Souter  Mean  Diameter  of  droplets  emerging 
from  the  nozzle  during  the  course  of  the  injection  has  been  evaluated  using 


equation  (3.1)  in  the  form 


,  -0 .1  35  ,  .0.121  .  ,0.131 

ds  =  A (AP )  (p  )  (Q) 

'a  m 


where , 


A  =  23.9  for  a  hole  nozzle 


2 


remain  in  the  core  of  the  spray,  the  smaller  ones  are  considered  to  adhere 
to  the  outer  fringes  of  the  spray. 

(3)  For  the  no-swirl  case,  the  spray  is  considered  to  be  a  right 
circular  cone  subtending  a  constant  angle  at  the  apex.  This  is  a  fairly 
safe  assumption  for  L/d  ratios  of  4  or  slightly  less.  Varde,  et  al  [15] 
have  shown  (Fig  2)  that  nozzles  with  L/d  ratios  of  4  do  not  exhibit  any 
change  in  spray  cone  angle  with  change  in  pressure  differential  across  the 
nozzle.  For  nozzles  with  L/d  ratio  less  than  4,  an  increase  in  gas  density 
or  AP  across  the  nozzle  tends  to  increase  the  spray  cone  angle. 

Therefore,  in  actual  engines  with  the  piston  moving  up,  the  effect  of  a 
decrease  in  pressure  differential  on  the  spray  cone  angle  is  nullified  by 
the  reverse  effect  due  to  a  simultaneous  increase  in  gas  density.  The 
value  of  the  constant  spray  angle  has  been  computed  using  the  expression 
given  by  Reitz  and  Bracco  [32] . 

Pa 

Tan  0/2  =  0.1  3  (1  +  — )  (4.3) 

Pi 

where , 

Pa  =  Chamber  gas  density 
=  Liquid  fuel  density 

For  purposes  of  convenience  in  computations,  and  to  facilitate  the 
study  of  spatial  and  temporal  distribution  of  liquid  and  vaporized  fuel 
within  the  spray,  the  spray  cone  has  been  subdivided  into  4  concentric 
subcones  under  conditions  of  no  swirl.  The  innercone  subtends  an  angle 
of  0/4  at  the  apex  and  is  a  solid  cone  while  the  second,  third  and 
fourth  cones  subtend  angles  2  0/4,  3  0/4  and  0  at  the  apex  respectively 


and  are  hollow  cones  which  envelope  the  inner  cones. 


(4)  The  distribution  of  liquid  fuel  by  moss  in  each  of  the  subcones 
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is  assumed  to  be  Gaussian.  As  we  move  away  from  the  inner  cone  the  amount 
of  fuel  contained  in  each  of  the  subcones  decreases. 


(5)  The  evaporation  rate  constant,  Kq  ,  for  diesel  fuel  in  stagnant 


air  was  obtained  by  interpolating  the  experimental  values  given  in 


Reference  5.  The  rate  at  which  liquid  fuel  evaporates  is  assumed  to  follow 

2 


the  "D  law”  which  relates  the  final  size  of  a  droplet  to  the  initial  size 
by: 

2  2 

(final  diameter)  =  (initial  diameter)  -  Kt  (4.4) 

where , 

K  is  the  evaporation  rate  constant  and  t  is  the  time  in  which  the 
change  in  diameter  takes  place. 


Here,  the  evaporation  rate,  K,  is  different  from  K  mentioned  earlier 

o 


as  the  ambient  gas  in  the  engine  cylinder  is  not  stagnant.  The  two 
constants  are  related  by  the  following  expression 


K  =  K  (1  +  0.276  Re^  Sc°’33) 
o 


(4.5) 

where , 

Re  =  Reynolds  number 

Sc  =  Schmidt  number,  taken  to  be  unity 

(6)  Since  there  is  relative  motion  between  droplets  and  the 
surrounding  air,  the  droplets  experience  a  drag  force  which  is  computed 


where , 


C  =  Coefficient  of  drag 
D  y 

d  =  Mass  averaged  diameter  of  droplet 
p  =  Ambient  gas  density 

d 

V  *  Relative  velocity  between  fuel  droplets  and 
surrounding  air  particles 

The  drag  force  is  composed  of  viscous  drag  and  pressure  drag  which 
depends  primarily  on  the  Reynolds  number  and  the  rate  of  mass  efflux  from 
the  droplet  surface.  The  coefficient  of  drag  to  be  used  in  equation  (4.6) 
can  be  approximated  as  a  function  of  the  Reynolds  number  and  is  given  by 


CD  =  (24/Re)  (1  +  0.15  Re0,687) 


(4.7) 


Yuen  and  Chen  [31],  based  on  their  experimental  results,  suggested 
use  of  the  above  equation  for  the  coefficient  of  drag  provided  the  Reynolds 
number  is  evaluated  using  the  density  of  the  free  stream  and  the 
characteristic  viscosity  coefficient  is  calculated  at  a  reference 
temperature  given  by: 


where , 


T  =  T,  +  (T  -  T,  ) / 3 
r  1  a  1 


=  Ambient  gas  temperature 


=  Fuel  Temperature 


(7)  The  initial  velocity  of  liquid  droplets  can  be  calculated  by 
applying  Bernoulli’s  equation  at  the  nozzle  exit  for  constant  injection 
pressure.  This  initial  velocity  is  assumed  to  be  the  same  for  all  droplets 


irrespective  of  their  size  and  may  be  found  from 

4 

V.  =  C .  (2g  (p .  .  -  p  . )/p,  ) 

lo  d  in]  cyl  1 


(4.8) 


where , 


Coefficient  of  discharge  for  the  nozzle 


P  =  Injection  pressure 

in  j 


P  ,  =  Cylinder  pressure 

cyl 


p^  =  Liquid  fuel  density 


(8)  The  gravity,  bouyancy,  and  pressure  forces  are  very  small 


compared  to  the  acceleration  and  drag  forces  and  as  such  have  been 


neglected.  For  the  same  reason,  Newton's  second  law  of  motion  in  its 


simple  form  can  be  used  to  trace  the  history  of  motion  of  fuel  droplets 


d/dt  (m.  •  V.  )  =  1/2  C.  A.  p  V 

1  1  D  1  a  r 


(4.9) 


V  =  V  expi-hp  A. C  )  (X) 
r  or  Ka  1  D 


(4.10) 


where , 


Mass  of  liquid  droplet 


V  =  Velocity  of  liquid  droplet 


Coefficient  of  drag 


Surface  area  of  droplet 


■  Ambient  gas  density 


V  *  Relative  velocity  of  droplet 


V  =  Initial  relative  velocity  of  droplet 
or 


(9)  Air  swirl  has  been  assumed  to  be  a  solid  body  rotation  around 


the  cylinder  axis  [8]  whereas  the  air  motion  due  to  piston  movement  and 


squish  has  been  neglected. 


mm 


(10)  Pressure  equilibrium  is  maintained  throughout  the  cylinder, 
i.e.  changes  in  the  engine  cylinder  pressure  are  assumed  to  stabilize 
instantaneously . 

(11)  Collisions  among  droplets  and  the  impingement  of  droplets  on 
the  combustion  chamber  wall  have  not  been  analyzed. 

(12)  Swirl  causes  the  concentric  right  circular  cones  to  be  deformed 
into  non-concentr ic  cones  with  elliptical  bases.  The  minor  axes  of  these 
cones  being  equal  to  the  diameter  of  the  cones  without  swirl. 

(13)  Vapors  generated  from  droplets  in  motion  form  a  homogeneous 
mixture  with  the  surrounding  air  and  move  with  the  same  angular  and  linear 
velocity  as  the  surrounding  air. 

(14)  The  rapid  pressure  rise  period  follows  immediately  after  the 
end  of  ignition  delay  which  is  given  by  equation  (3.2).  The  values  of  the 
constants  have  been  taken  from  reference  5  to  give  equation  (3.2)  the  form 

t  =  0.872  (P  )_1  *24($)-2'1  °exp(4050/T  )  (4.11) 

a 

where , 

P^  =  Pressure  of  gas  in  chomber  in  atm 

<J>  =  Oxygen  concentration  in  ambient  gas 

T  «  Mean  temperature  of  ambient  gas  in  °K 
a 

(15)  Combustion  begins  in  a  zone  after  the  end  of  ignition  delay 
only  if  the  equivalence  ratio  lies  within  the  flammability  limits  of  0.3  to 
3.0. 

(16)  All  the  premixed  fuel  identified  to  ignite  burns  in  three 
degrees  of  crank  angle  rotation  [23], 


07) 


After  the  end  of  the  rapid  pressure  rise  period,  the  fuel  is 


assumed  to  burn  at  the  rate  it  is  injected  with  equivalence  ratio  equal  to 


unity . 


4.2  MODEL  INPUTS 

The  following  are  the  necessary  inputs  to  the  model: 

(1)  Spray  cone  angle 

(2)  Sauter  Mean  Diameter  of  fuel  droplets  and  their  distribution 

(3)  Mass  averaged  diameter  of  droplets 

(4)  Liquid  fuel  mass  distribution  across  the  spray 

(5)  Fuel  properties  such  as  evaporation  rate  constant,  density, 
and  boiling  point  temperature  as  functions  of  pressure  and 
temperature 

(6)  Fuel  injection  pressure 

(7)  Swirl  ratio 

(8)  Engine  geometry 

(9)  Engine  speed 

(10)  Fuel  injection  timing 

(11)  Duration  of  injection 

(12)  Injector  nozzle  diameter  and  its  coefficient  of  discharge. 


4.3  SUMMARY  OF  MODEL 

The  model  requires  fuel  properties  and  the  engine  dimensions  as  the 
main  input  parameters  and  analyzes  the  spray  as  an  agglomerate  of  tiny 
droplets,  forming  a  heterogeneous  mixture  and  existing  under  supercr itical 
conditions . 
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A  finite  difference  numerical  approach  has  been  adopted  that  creates 
new  spray  segments  with  the  passage  of  time  between  the  stations  traversed 
by  the  fuel  droplets  in  small  but  equal  intervals  of  time.  The  sizes  of 
segments  vary  with  time  and  also  differ  from  cone  to  cone  depending  on  the 
distance  travelled  by  different  sized  droplets  in  the  same  time  interval. 
The  initial  momentum  of  the  spray  in  each  cone  at  the  start  of  injection  is 
the  momentum  of  the  fuel  injected  in  the  first  time  increment.  As  the 
droplets  move  through  the  air  they  experience  drag  and  evaporation  and  at 
the  end  of  a  time  interval  the  total  momentum  is  shared  by  the  remaining 
liquid  fuel,  vaporized  fuel  and  the  entrained  air.  An  identical  pattern  is 
followed  oy  the  remaining  liquid  droplets  in  the  next  increment  of  time  and 
also  by  the  succeeding  droplets  coming  out  in  the  following  interval  of 
time.  Thus,  the  basic  equations  of  mass  and  momentum  are  repeatedly 
satisfied  to  solve  for  decaying  velocity  of  droplets  between  stations  and 
the  increasing  velocity  of  air  in  the  direction  of  penetration.  The  non- 
Isothermal  character  of  the  spray  is  recognized  by  simultaneously 
establishing  an  energy  balance  in  each  time  step.  Such  a  procedure  keeps 
track  of  droplet  history,  the  amount  of  fuel  vaporized  and  the  quantity  of 
air  entrained  in  a  particular  time  step  at  different  locations  of  the 
spray.  This  enables  the  model  to  specify  the  temporal  and  spatial 
variation  of  combustible  mixture  in  various  regions  of  space. 

Another  salient  feature  of  the  model  is  its  capacity  to  take  into 
account  the  swirling  of  air  and  the  variations  in  back  pressure,  chamber 
volume  and  chamber  gas  temperature  with  time  due  to  the  continuous  motion 
of  tie  piston  in  the  cylinder.  The  progressive  deformation  of  the 
concentric  cones  because  of  swirl  has  been  followed  by  determining  the 


deformed  boundaries  defined  by  the  locus  of  the  positions  of  droplets  in 


each  cone  at  the  end  of  each  time  interval  and  the  corresponding  deflection 
of  fuel  vapors  with  swirling  air  in  the  same  interval  of  time. 

Figure  4  shows  the  velocity  and  displacement  vector  diagram  for  a  fuel 
droplet  in  the  first  three  increments  of  time.  1-2’  is  the  intended  line 
of  motion  of  the  droplet  in  the  first  increment  of  time.  Point  2  is  the 
actual  location  of  the  droplet  at  the  end  of  the  first  time  interval  due  to 
superposition  of  the  swirl  velocity.  The  figure  also  shows  the  location  of 
the  droplet  at  points  3  and  4  in  the  next  two  time  intervals.  The 
penetration  of  the  droplet  at  the  end  of  the  "n  th  "  increment  of  time  is 
given  by  PEN  (n)  and  is  different  from  the  distance  actually  travelled  by 
the  droplet  along  a  curved  path.  Vapors  generated  are  assumed  to  have  the 
same  linear  and  angular  velocity  as  that  of  the  surrounding  air  which  makes 
the  fuel-air  distribution  at  the  end  of  the  3rd  interval  of  time  appear  as 
shown  in  Figure  5.  The  final  configuration  of  the  spray  is  depicted  in 
Figure  6  with  elongation  of  the  circular  base  taking  place  in  the  direction 
of  swirl.  Since  there  is  no  force  considered  to  be  acting  perpendicular  to 
the  plane  of  the  spray,  no  deformation  takes  place  in  that  direction.  This 
is  primarily  the  reason  why  the  minor  axes  of  the  elliptical  bases  of  the 
deformed  cones  are  considered  to  be  equal  to  the  diameter  of  the  base  of 
the  undeformed  cone  without  swirl. 

Keeping  account  of  the  spray  geometry  at  all  times  makes  it  possible 
to  detemine  the  intermixing  of  fuel  vapors  and  air  between  overlapping 
cones  when  computing  the  equivalence  ratio  at  various  locations  within  the 
fuel  spray. 
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Fig.  4.  Velocity  and  Displacement  Diagram 


Figure  7  gives  an  idea  of  how  the  movement  of  vapors  takes  place  from 
outer  to  inner  cones.  While  in  the  case  of  no  swirl,  the  vapors  generated 
in  a  cone  may  be  assumed  to  be  confined  to  the  annular  region  between  the 
cones,  the  same  assumption  cannot  be  extended  to  the  case  with  swirl. 

Since  the  penetration  of  droplets  in  the  outer  cone  is  not  the  same  as  that 
of  inner  cones  and  the  size  of  the  segments  varies  also,  one  or  more 
segments  of  outer  cones  may  contribute  vapors  to  the  inner  cones.  Also, 
since  one  or  more  of  these  contributing  segments  may  not  be  formed  at  the 
time  of  observation,  the  change  in  equivalence  ratio  in  any  segment  due  to 
overlapping  varies,  making  equivalence  ratio  at  a  particular  location  a 
function  of  time.  The  mass  transfer  of  vapor  from  segments  of  outer  cone 
to  that  of  inner  cones  is  considered  proportional  to  the  percentage  of 
overlapping  by  volume. 

Ignition  delay  for  each  time  dependent  segment  is  calculated 
individually.  Once  the  ignition  delay  is  over,  all  segments  between  the 
flammability  limits  of  0.3  to  3.0  burn  during  a  period  of  3  degrees  of 
crank  angle  rotation.  The  remaining  fuel  burns  at  the  same  rate  at  which 
it  is  injected  with  an  equivalence  ratio  of  unity.  A  program  developed  at 
the  University  of  Wisconsin  [29]  for  calculating  properties  of  products  of 
combustion  is  called  as  a  subroutine  to  determine  the  heat  released  by  each 
segment  and  the  resulting  adiabatic  flame  temperature  attained  by  that 
segment.  Pressure  rise  is  finally  determined  by  equating  the  product  of 
instantaneous  cylinder  pressure  and  volume  to  the  summation  of  mRT  of 
burned  segments,  unburned  segments  and  surrounding  air  at  the  end  of  each 
time  step.  A  listing  of  the  computer  program  is  given  in  Appendix  1. 
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5.  DISCUSSION  OF  RESULTS 


5.1  COMPARISON  OF  MODEL  RESULTS  WITH  EXISTING  DATA 

The  prediction  of  the  model  for  spray  penetration  with  zero  swirl 
ratio  has  been  compared  in  Figures  8  and  9  with  data  obtained  by  Takeuchi, 
et  ol  [16]  and  Varde  and  Popa  [30]  in  two  seperate  experiments  conducted  in 


1 


quiescent  air.  Every  care  has  been  taken  to  extract  information  concerning 
inputs  to  the  model  from  References  16  and  30  in  order  to  make  the 
comparison  justifiable.  Essential  information  on  inputs  known  to  have 
substantial  bearing  on  the  outcome  of  results  if  not  available  in  the  above 
references  have  had  to  be  assumed.  Lack  of  knowledge  on  either  one  or  more 
of  the  following: 

( 1  )  Nozzle  coefficient  of  discharge 

(2)  Nozzle  L/d  ratio 

(3)  Spray  cone  angle 

(4)  Properties  of  fuel  used 

(5)  Chamber  gas  properties 

(6)  Breakup  time  or  breakup  length  of  liquid  fuel  element 
appears  to  be  the  major  factor  responsible  for  small  discrepancies  noticed. 

Not  much  published  data  on  penetration  observed  on  injecting  fuel  into 
air  in  motion  is  available.  Hiroyasu,  et  al  [13]  have  recommended  the  use 
of  a  factor  which  is  a  function  of  swirl  ratio  to  evaluate  the  value  of 
penetration  in  swirling  air  provided  the  value  of  penetration  in  stagnant 
air  is  known.  This  factor  has  been  applied  to  Hiroyosu’s  experimental  data 
available  on  penetration  without  swirl  [13]  and  plotted  against  the  models 
prediction  in  Figure  10.  Figure  11  shows  a  similar  comparison  between 
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COMPARISON  WITH  HIROYASU'S  MODEL  WITH  SWIRL 


timo  ms 

COMPARISON  WITH  UNIVERSITY  OF  MANCHESTER  MODEL  IN  CROSS  FLOW 


calculated  values  and  the  experimental  values  of  Yule,  et  al  [19]  obtained 

in  a  wind  tunnel  under  constant  back  pressure.  The  calculated  penetration 

is  slightly  higher  to  begin  with  but  becomes  less  than  the  experimental 

value  1  ms  from  the  start  of  injection.  It  may  be  pointed  out  again  that 

spray  penetration  with  swirl  as  given  by  the  model  is  the  distance  of  the 

last  location  of  the  leading  droplet  from  the  nozzle  tip.  It  is  different 

from  the  actual  distance  travelled  by  the  droplet  along  a  curved  path. 

5.2  EFFECT  OF  SOME  ENGINE  OPERATING  VARIABLES  ON  PENETRATION  AND 
EQUIVALENCE  RATIO 

5.2.1  EFFECT  OF  INJECTION  PRESSURE 

Figures  12  and  13  show  the  calculated  effect  of  injection  pressure  on 
spray  penetration  for  constant  load.  Figure  12  has  been  plotted  from 
results  obtained  by  keeping  the  back  pressure  constant  at  a  value 
corresponding  to  engine  cylinder  pressure  at  the  start  of  injection,  viz. 

24  atm.  Figure  13  shows  a  similar  variation  under  engine  running  con¬ 
ditions.  The  comparison  between  the  two  conditions  for  injection  pressures 
of  220  atm  and  100  atm  has  been  made  in  Figure  14.  The  penetration  in  both 
cases  increases  with  an  increase  in  injection  pressure  but  the  penetration 
under  engine  conditions  is  less  at  the  same  injection  pressure  due  to  the 
increased  drag  imposed  by  the  denser  air.  The  constant  back  pressure  case 
has  been  included  for  comparison  purposes  since  most  of  the  experimental 
data  available  was  obtained  in  constant  pressure  chambers. 


TIME  ma 

FIG.  12  EFFECT  OF  INJECTION  PRESSURE  ON  SPRAY  TIP  PENETRATION 


TIME  ms 

FIG.  14  COMPARISON  BETWEEN  PENETRATION  UNDER  CONSTANT  BACK  PRESSURE  AND  ACTUAL  ENGINE  CONDITION 


The  distribution  of  air  fuel  mixture  in  various  subcones  after  about 
1  ms  from  the  start  of  injection  for  both  the  running  engine  condition  and 
the  constant  back  pressure  condition  are  shown  using  bargraphs  in  Figures 
15  and  16.  The  time  1  ms  corresponds  approximately  to  the  end  of  ignition 
delay  under  typical  engine  operating  conditions.  The  height  of  the  bar  at 
any  time  "t"  describes  the  equivalence  ratio  in  a  segment  of  a  cone  formed 
at  the  end  of  ignition  delay  from  vaporizing  droplets  that  were  injected 
"t"  ms  before  the  end  of  ignition  delay.  Thus,  for  example,  a  group  of 
bars  at  0.22  ms  represent  equivalence  ratio  in  cone  segments  close  to  the 
nozzle  tip  at  the  end  of  ignition  delay  while  those  at  time  0.98  ms  repre¬ 
sent  cone  segments  near  the  spray  tip  at  the  same  instant.  With  an  in¬ 
crease  in  injection  pressure  at  constant  engine  load,  the  duration  of 
injection  decreases  and  a  greater  mass  of  fuel  is  injected  in  the  same  time 
interval.  However,  this  does  not  result  in  higher  values  of  equivalence 
ratio.  Indeed,  the  value  decreases  due  to  the  fact  that  penetration  is 
greater  at  higher  injection  pressures  which  allows  for  more  air 
entrainment.  Also,  due  to  higher  temperatures  and  smaller  penetrations  one 
would  expect  richer  fuel  air  mixtures  in  the  case  of  running  engine  condi¬ 
tions  as  compared  to  constant  back  pressure  conditions.  This  trend  is 
noticed,  but  only  near  the  nozzle  tip,  (Figure  17).  Rapid  evaporation  of 
droplets  in  the  initial  stages  leaves  less  liquid  fuel  farther  away  from 
the  nozzle  which  explains  the  occurrence  of  this  effect.  It  is  interesting 
to  note  that  ignition  delay  is  not  altered  by  changing  the  injection 
pressure  because  of  the  weak  influence  of  equivalence  ratio  on  ignition 


delay  period. 


UNDE*  ACTUAL  ENGINE  CONDITIONS 
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5.2.2  EFFECT  OF  INJECTION  TIMING 

Spray  tip  penetration  decreases  as  injection  timing  is  retarded  due  to 
the  increase  in  chamber  gas  density  as  shown  in  Figures  18  and  19.  The 
reduction  in  penetration  is  more  in  the  case  of  actual  engine  conditions  as 
compared  to  constant  back  pressure  conditions  for  the  same  change  in  injec¬ 
tion  timing  as  shown  in  Figure  20. 

The  variation  of  equivalence  ratio  with  injection  timing  is  shown  in 
Figures  21  and  22.  For  a  constant  duration  of  injection,  the  distribution 
of  fuel  air  mixture  is  slightly  leaner  in  the  case  of  retarded  injection 
timing.  This  is  probably  due  to  higher  back  pressures  associated  with 
retarded  injection  which  results  in  smaller  droplet  velocity  and  a  smaller 
quantity  of  fuel  injection  in  the  same  interval  of  time,  inspite  of  the 


fact  that  penetration  and  air  entrainment  is  also  reduced.  In  addition, 
in  the  case  of  engine  operating  conditions,  the  ignition  delay  period  is 
also  reduced  with  retarded  injection  thereby  reducing  the  quantity  of 
total  fuel  injected  during  the  ignition  delay  period.  This  explains  why 
a  comparatively  leaner  distribution  is  noticed  in  the  case  of  actual 
engine  conditions. 
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EFFECT  OF  INJECTION  TIMING  ON  SPRAT  TIP  PENETRATION 
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EFFECT  OF  INJECTION  TIMING  ON  EQUIVALENCE  RATIO  DISTRIBUTION  IN  DIESEL  SPRAYS 
UNDER  CONDITIONS  OF  SWIRL  WITH  BACK  PRESSURE  HELD  CONSTANT 
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5.2.3  EFFECT  OF  ENGINE  SPEED 

There  is  not  much  change  noticed  in  the  spray  tip  penetration  as 
the  engine  speed  is  varied  provided  the  total  load  on  the  engine  is 
maintained  constant.  This  was  achieved  in  the  model  by  increasing  the 
total  duration  of  fuel  injection  and  thus  allowing  for  the  same  quantity 
of  fuel  to  be  injected.  The  nominal  effect  of  engine  speed  on  penetration 
is  shown  in  Figure  23. 

5.3  PRESSURE  RISE  DIAGRAM 

As  mentioned  earlier,  the  pressure  rise  following  ignition  delay  has 
been  evaluated  on  the  assumption  that  the  "spike"  lasts  for  3  degrees  of 
crank  angle  and  thereafter,  the  fuel  burns  at  the  rate  it  is  injected  with 
equivalence  ratio  of  1 .  A  cylinder  pressure  versus  crank  angle  diagram  is 
shown  in  Figure  24.  The  curve  a-b-c-d  describes  the  pressure  variation 
during  motoring.  The  curve  a-b-c ' -d ’ -e-f-g  represents  the  pressure 
changes  in  the  cylinder  when  combustion  takes  place.  In  this  second 
curve,  b-c’  corresponds  to  the  ignition  delay  period,  the  premixed  fuel 
burns  between  c'  -d',  the  remaining  fuel  is  consumed  by  diffusion  flame 
between  d'  -e  and  finally  e-f-g  is  the  outcome  of  compression  and 
expansion  due  to  piston  movement. 
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6.  EXPERIMENTAL  SETUP 


6.1  ENGINE  DESCRIPTION 

A  Ricardo  Hydra  Series  Single  Cylinder  Direct  Injection  Research 
Diesel  Engine  has  been  installed  in  the  Engine  Laboratory.  Experiments  have 
been  done  on  this  engine  to  validate  the  results  given  by  the  model.  The 


term  "running  engine  conditions"  in  the 
engine  in  particular . 

Some  of  the  salient  features  of  th 
Engine  details: 

Bore  size 
Stroke  length 
Rated  speed 

Peak  cylinder  pressure 
Measured  compression  ratio 

By  changing  the  crankshaft 
varied  between  71  mm  to  105  mm. 
Fuel  injection  system: 

Nozzle 

Nozzle  Opening  Pressure  (NOP) 
Instrumentation : 

Engine  speed 
Engine  load 
Fuel  flow 
Air  flow 

Nozzle  needle  lift 


preceding  chapters  refers  to  this 

s  engine  are: 

80 . 26  mm 
88.90  mm 
75  rev/sec 
120  bars 
20/1 

I  connecting  rod,  the  stroke  can  be 

4  holes,  0.21  mm  diameter 
250  bars 

digital  display 

load  cell  with  digital  display 
AVL  gravimetric  flow  meter 
Miriam  Laminar  flow  element 
inductive  transducer 


Fuel  line  pressure 


strain  guage  transducer 


Cylinder  pressure 


Kistler  piezoelectric  transducer 


Air  inlet  temperature 
Water  outlet  temperature 
Oil  inlet  temperature 
Fuel  temperature 
Exhaust  outlet  temperature 
Valve  timing  : 

Inlet  opens 
Inlet  closes 
Exhaust  opens 
Exhaust  closes 


thermocouple 

thermocouple 

thermocouple 

thermocouple 

thermocouple 

10  BTDC 
42  ABDC 
58  BBDC 
10  ATDC 


Recommended  static  injection 
Speed  rev  /  sec 
20 
30 
40 
50 
60 


timing  (S.I.T.)  under  full  load 
S.I.T.  BTDC 
13 
16 
17 
19 
22 


Restrictions : 

Maximum  exhaust  temperature 

Lube  oil  inlet  temperature 

Water  outlet  temperature 

Exhaust  back  pressure 

Inlet  air  temperature 

Inlet  fuel  temperature 
(to  pump) 


750  C 
85  C 
85  C 
0 . 1  bar 
18-22  C 
15-20  C 
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6.1.1  AIR  INLET  SYSTEM 

The  inlet  system  consists  of  a  cast  aluminium  alloy  manifold  provided 
with  a  1  KW  air  heater,  a  Miriam  50  MC2  laminar  flow  element  and  an  air 
filter.  The  temperature  of  air  passing  through  the  heater  can  be 
controlled  by  an  air  temperature  controller  at  the  control  console.  The 
signals  for  this  controller  are  provided  by  a  duplex  thermocouple. 

6.1.2  FUEL  SYSTEM 

The  system  consists  of  a  lift  pump  which  is  supplied  fuel  from  the 
fuel  reservoir  and  delivers  it  to  a  Bosch  injection  pump  for  injection  into 
the  cylinder.  The  injector  is  detachable  and  can  be  changed  to  study  the 
effect  of  injector  dimensions  on  engine  performance .  There  is  a 
provision  for  recirculating  excess  fuel  from  the  injector  pump  back  to  the 
flowmeter.  The  injector  unit  is  provided  with  a  thermocouple  to 
measure  the  inlet  fuel  temperature  which  can  be  monitored  from  the  control 
console.  Fuel  injection  pressure  is  monitored  with  a  strain  gauge  type 
pressure  transducer.  An  Inductive  Displacement  transducer  is  used  to 
measure  the  needle  lift.  Its  signal  is  conditioned  by  a  carrier  amplifier 
for  display  on  an  oscilloscope.  The  quantity  of  fuel  delivered  per  cycle 
can  be  remotely  controlled  by  adjusting  the  throttle  and  the  injection 
timing  can  be  set  by  remote  control  between  0  and  40  BTDC. 

An  AVL  730  gravimetric  fuel  consumption  measuring  system  is  also 
fitted  in  the  fuel  circuit.  The  change  in  weight  of  the  vessel  containing 
fuel  is  determined  by  the  change  in  its  position  over  a  frictionless  blade 
spring.  The  change  is  detected  by  a  highly  sensitive  capacitive  displace¬ 
ment  transducer.  The  equipment  offers  the  possibility  to  compute  and  indi- 
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(a)  cumulative  consumption  over  a  preset  measuring  period 

(b)  Mean  consumption  rate  over  a  selected  measuring  period 

(c)  The  instantaneous  mass  flow  rate. 

6.1.3  COOLING  SYSTEM 

The  closed  circuit  pressurized  engine  coolant  system  uses  water  and 
antifreeze.  It  is  used  to  cool  both  the  engine  block  and  the  cylinder 
head.  A  duplex  thermocouple  at  the  header  tank  controls  the  flow  rate  of 
cooling  water.  The  cooling  water  is  also  passed  through  a  heat  exchanger 
for  cooling  the  lubricating  oil. 

6.1.4  OIL  SYSTEM 

The  lubricating  system  is  used  both  as  a  lubricating  unit  and  a 
cooling  unit.  Oil  under  pressure  is  supplied  to  crankshaft  bearings,  big 
end  bearings,  and  camshaft  bearings.  A  safety  system  gives  a  warning 
signal,  upon  loss  of  oil  pressure  or  if  oil  temperature  increase  beyond 
85  C. 

6.1.5  PRESSURE  MEASUREMENTS  AND  DATA  ACQUISITION 

Cylinder  pressure  is  monitored  with  a  Kistler  Piezoelectric 
transducer.  This  transcucer  transmits  signals  through  a  charge  amplifier 
to  the  oscilloscopes.  A  Data  Precision  Data  6000  4  channel  digital 
oscilloscope  is  used  to  monitor  cylinder  pressure,  needle  lift  and 
crankshaft  position.  The  Data  6000  4  channel  plug  in  has  a  14  bit  100  KHz 
A/D  converter  (25  KHz  per  channel)  and  can  process  and  store  the  data  on 
floppy  disks  or  interface  to  a  computer. 


7.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  outcome  of  this  study  has  been  the  development  of  a  mathematical 
model  to  represent  spray  development,  heat  release  and  pressure  time 
history  in  direct  injection  diesel  engines.  The  model  is  capable  of 
predicting  : 

(i)  Spray  tip  penetration 

(ii)  The  local  fuel  air  distribution  within  the  spray 

(iii)  The  ignition  delay  period 

(iv)  Locations  in  spray  where  autoignition  begins 

(v)  The  pressure  rise  due  to  fuel  combustion 

under  continuously  changing  environmental  conditions  as  may  be  encountered 
in  actual  engines. 

In  addition  to  predicting  the  above,  the  model  can  effectively  be 
employed  to  model  the  effects  of  varying  some  of  the  engine  operating 
variables  such  as  injection  pressure,  injection  timing,  engine  speed, 
nozzle  dimensions,  etc.  on  each  of  the  above  characteristics  of  a  diesel 
spray . 

The  commonly  available  inherent  characteristics  of  an  engine  and  the 
properties  of  the  fuel  serve  as  the  bulk  of  the  essential  inputs  required 
by  the  model.  An  additional  feature  of  the  model  is  its  capability  to 
determine  the  geometry  of  the  spray  under  conditions  of  swirl. 

Although  experiments  have  yet  to  be  conducted  to  validate  the  efforts 
of  this  exercise,  the  results  of  the  model  hove  been  compared  to  the 
existing  data  from  experiments  performed  in  constant,  volume  hot  and  cold 
bombs  and  found  to  be  in  good  agreement. 


cylinder  pressure  history  with  engine  experiments  using  hexadecane  fuel, 


they  will  be  forwarded  as  an  addendum  to  this  report. 
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APPENDIX  1 

0** ********************************************* ****************** 

C  THIS  PROGRAM  PREDICTS  THE  SPRAY-TIP  PENETRATION  AND  EQUIVALENCE 
C  RATIO  DISTRIBUTION  WITHIN  THE  GEOMETRY  OF  SPRAY  IN  DIRECT 
C  INJECTION  DIESEL  ENGINE  UNDER  CONDITIONS  OF  SWIRL. 

0***** ************************************************** ********** 

REAL  VX(61,61) ,SI(60) ,XTT(60) ,ATHXV(60)  ,XTPV(60) ,FTHXV(60) 

REAL  VL1 (61,61) ,VR1 (61,61) ,XRE1 (61,61) ,CD1(61,61) 

REAL  XMAIR(4,61) ,FFR(61)  ,XXTP(6I) ,YXTP(61) ,XXTPV(61) ,YXTPV(61) 

REAL  XK0{5) ,XK2 (61,61) ,XKK(61) , PEN (4, 61) ,XC1(61,61) ,XK1(61,61) 

REAL  FR(5) ,DSS(5) ,DM(5) ,TH(5) 

REAL  TFL{61,61) ,TVOLL (61, 61) ,CV(61) ,CL(61) 

REAL  BCD (6 1,61) 

REAL  YXT (61) ,XXFV(61) ,SUMR(61) ,TSUM(61) ,FMASS(61) 

REAL  D(61, 61) ,AVGD2 (61, 61) ,XM1(61) ,XND(61) ,VR2(61,61) 

REAL  DS (61,61) ,AVGDS (61,61) ,XNDS(61) ,XMS1(61) ,AVGVR2(61,61) 

REAL  XT (61) ,XRE2 (61,61) ,CD2(61,61) ,AVGCD2 (61,61) ,XC2 (61,61) 

REAL  FV(61,61) ,FL(61,61) ,XMA(61,61) ,AVGXC(61,61) 

REAL  VR(61)  ,VOLFV (61,61) ,VOLL(61,61) ,VOLA(61,61) 

REAL  AA(0:61) ,AM(61,61,61)  ,PHE(4,61) 

REAL  V ( 61)  ,X(61,61) 

REAL  Xl(61,61) 

REAL  TIM (60) , DELAY (4 , 61 , 61) 

REAL  FFW(61)  ,XXMA(61)  ,DUMMY1(61)  ,VOA(61) 

REAL  AVOA (61) ,FAR(4,61,61) ,TFV(4, 0:61, 0:61) ,TA(4,61,61) 

REAL  AVOAI ( 10) ,FTA (4 , 61, 61) ,TAC (4, 61, 61) ,TAM(4,61,61) 

COMMON/ A/ AREA ( 2) 

DOUBLE  PRECISION  THXL , ARXL ( 60 ) , XLAMDA1 , XLAMDA2 , XXLAMDA , BETA1 , BETA2 
DOUBLE  PRECISION  BETA, SIGMA, VIJ ,BTA, SV,AI , RXL,XTP (61) ,DUMMY2 (61) 
DOUBLE  PRECISION  ARR 

(2***************************** ***★★*★****★**★ ******** **★★*★****★★★**** 

C  INPUT  DATA  :  THIS  NEEDS  TO  BE  CHANGED  WITH  EVERY  TEST  RUN  TILL  SUCH 
C  TIME  THAT  WE  OBTAIN  A  GENERAL  EXPRESSION  FOR  EACH  ONE  OF  THEM. 

Qkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk  kirk  it  It  kkkkkkk  ★★★***★* 

DATA  FR/ 0., .60, .25, .10,. 05/ 

DATA  XK0/0. , .009, .011, .012, .015/ 

DATA  DSS/0 . , 80 . E-4 , 56 . E-4 , 44 . E-4 , 20 . E-4/ 

DATA  DM/0.,48. E-4, 36. E-4, 24. E-4, 12. E-4/ 

DATA  TH/0.,.026,. 052,. 078,. 105/ 

Qk  k  kirk  kirk  irk  it  kirk  irk  it  kirk  it  it  it  kit  Irk  it  irk  irk  kirk  it  it  it  it  kirk  k  kirk  kk  it  it  it  it  irk  it  it  irk  kit  irk  kirk  kkkk 

AVG (X, Y) = (X+Y) /2 . 

FVL (PINJ, PA) =. 4*SQRT( 1962000.* ( PIN J- PA) /. 74871) 

FXRE (VR, D , CONST) =VR*D*CONST 
FCD(XRE) =24 ./XRE* ( 1 . + . 15*XRE** .687) 

FXK (XK0 , XRE, SC) =XK0* ( l.+. 276*SQRT (XRE)  *SC**0. 333334) 

FD (DOLD, XK , DT) = (DOLD*DOLD-XK*DT) 

FXC (DNEW, CD, CONST) =4 ./3 .*DNEW/CD*CONST 
FMFL(N,D) =3. 1416/6. *N*D**3. 

FVR(V,X1,X2,XC) =V*EXP (- (X2-X1) /XC) 

OPEN  (UNIT=1,FILE='  SPRAY.DAT'  , STATUS= ’ NEW'  ) 
0********************************************************************* 

C  CONSTANTS  FOR  FUEL  TO  BE  USED  IN  THE  EXPRESSION  FOR  IGNITION  DELAY. 

C  SHOULD  BE  CHANGED  FOR  DIFFERENT  FUELS. 
0********************************************************************* 
A=0. 872 


fL- _ 


v.v.v; 


AN— 1.24 

C=—2.1 

DD=4050. 

PHID=1. 

0*  ************************************************************* ******* 

C  OTHER  KNOWN  INPUT  PARAMETERS  SOME  OF  WHICH  MAY  NEED  ALTERATIONS  WITH 
C  EACH  TEST  RUN. 

0*  ******************************************************************** 
PI=3. 146 
CDD= .  4 
DNN=.03 
R=2.8326 
TA0=305. 

TBP=480. 

GC=981. 

TL=480 . 

PINJ=220 

THINJ=11 

RAJ=70. 

RPM=900 . 

PT0=1. 0 
PFL=9. 5 
Z=.72 

RHOL=. 74871 
SC=1. 

FAS=. 066667 

0* *************************************************************  ***** 

C  INITIALIZE  VALUES 

C***** ************************************************************** 

VAW=0. 

XVAW=0 . 

XXVAW=0.0 
XT (1) =0. 

XTF=0 . 

XTP(l) =0. 

SI(1)=0. 

0* Hr*Hr**HrHrHrHrHrHrHrHrHr****HrHrHrHrHr***HrHrHr*HrHr*HrHrHrHcHr*HrHrHr*******HrHrHrHrHrHr  Hr  ★**★**★*★* * 

RHOA0=PT0/(R*TA0) 

PA0=PT0 

AG=.5 

IM12=60 

DT=THINJ/(IM12*6.*RPM) 

0***Hr*  Hr  Hr  Hr*  Hr  Hr  Hr  Hr  Hr  Hr  Hr  ★**★******★**★★★★★★★**★*★**★*★***★★★★★*★****★*★***** 

C  THIS  DO  LOOP  IS  RESPONSIBLE  FOR  THE  CREATION  OF  NEW  TIME  DEPENDENT 
C  SEGMENTS 

0Hr  Hr********HrHr*Hr*HrHrHrHrHrHr***Hr*Hr*HrHr**Hr*HrHrHrHr*HrHrHrHr*HrHr*HrHrHr****HrHr*HrHrHr******** 

DO  5000  LK= 1 , 60 

0*  ★★**★★★★★★★★★★*★★★**★★***★★*★★*★★★★★★*★★*★★★★*★**★***★★★★★***★*★*★* 

C  THIS  DO  LOOP  IS  FOR  ANALYSIS  IN  EACH  SUBCONE  TAKEN  ONE  AT  A  TIME 
0*  *★*★★★****★★**★★*★*★********★***★***********★*★*★★*******★★******** 

DO  111  IV=2,5 

0*  ******************************************************************* 

C  INITIALIZING  FOR  EACH  SUBCONE 

0*  ******************************************************************* 

VOA(l) =0.0 


AREA  ( 1)  =0 . 0 
TIME=0. 

TMASS=0. 

FVTOT=0 . 

VAW=0 . 

XVAW=0 . 

XXVAW=0.0 

XTV=0. 

ATHXV (1) =0. 

TFV (IV, 0,0) =0. 

TFV (IV, 0, 1) =0 . 

TFV (IV, 0,2) =0. 

PEN ( IV, 1) =0.0 

C*** ******************* ********************** ************************ 

C  SOME  KNOWN  ENGINE  CHARACTERISTICS-  LIQUID  FUEL  TEMP., INLET  PRESSURE, 
C  CYLINDER  VOLUME,  AMBIEOT  AIR  TEMP. ,  AND  FUEL  INJECTION  TIMING 

£***********★**★**************★******★********★************** ★★****★★* 

TL=480. 

P=l. 

V0=246. 

T=300. 

ANG=0. 38 

CALL  PV ( ANG, V0 , P,V1, PI) 

T1=TEMP(T,P,P1) 

TR=TL+ . 3333334* (Tl-TL) 

TM= (Tl-TL) /ALOG (Tl/TL) 

XMEU=10.* (7.13E-6+3.8E~8*TR) 

XK00=XK0 ( IV) 

£★*★*****★★***  *★****★***★★ *****************************************  *** 

C  THIS  DO  LOOP  ANALYSIS  EACH  OF  THE  PREVIOUS  SEGMENTS  EVERY  TIME  A  NEW 
C  SEGMENT  IS  FORMED 

c*  ★*★★**★★★★★★*★★★*******■*★**★★★★★★*★**★★★★*★*★★*★★*★★**********★**★★* 

DO  110  1=1, LK 

RHOA= (Pl-PFL) / ( Z*R*T1) 

IJ1=2 
TL=480 . 

D(l,  I)  =DM(IV) 

AVGD2 (1, I) =D(1, I) 

DS  ( 1, 1)  =DSS  ( IV) 

AVGDS  (1,1)  =DS  (1,1) 

TIME=TIME+DT 

J=I+1 

VL1 (1, I) =FVL(PINJ, PI) 

XMINJ=PI/4 .*DNN*DNN*VL1 ( 1 , I) *RHOL*DT*FR ( IV) 

TMASS=TMASS+XMINJ 
FL(1, I) =XMINJ 

XM1(I)=PI/6.*D(1,I)  **3.*RHOL 
XMS1 ( I) =PI/6 .*DS ( 1, 1)  **3.*RHOL 
XNDS (I) =XMINJ/XMS1 (I) 

XND(I) =XMINJ/XM1 (I) 

CONSTl=RHOA/XMEU 
CO  NST2=RHOL/RHOA 
IF(I.GT.l) GO  TO  141 
XT (J) =VLl (1,1) *DT+XT (J-l) 


GO  TO  142 

141  XT ( J) =VL1 ( I , 1-1) *DT+XT ( I) 

142  CONTINUE 

XTT ( J) =XT ( J) *COSD (RAJ) 

ANG=ANG-(RPM*2.* PI/60.) *DT 
CALL  PV (ANG, VI, PI, V2, P2) 

SM1=0. 

SM2=0. 

Q************************** ****** ********************************** 

C  INITIAL  ESTIMATE  OF  MASS  OF  AIR  TAKING  PART  IN  MOMENTUM  EXCHANGE 
III=IV-1 

HAR=PI/3.* (TAN (TH ( IV) ) *TAN(TH(IV) ) -TAN (TH (III) ) *TAN(TH(III) ) ) 
SMI = SMI + HAR* XT  (J) **3 
VOLA(J, I) =SM1 
XMA( J, I) =RHOA* SMI 
AIR=SM1*  RHOA 
11=1-1 

IF ( I . EQ. 1) GO  TO  52 
XMA(J, I) =XMA(J,I)-XMA(I, II) 

52  CONTINUE 

XMAIR ( IV, I) =XMA( J, I) 

N=I 

VR1  ( 1  / 1 )  =VL  1(1, 1 )  -XV AW 

0* *★***★*★★****★★***★**★*****★*★★********★**★★★*★***★*★***★★*★**★** 

C  THIS  DO  LOOP  REPEATS  THE  ANALYSIS  OF  NEXT  DO  LOOP  TO  APPROACH  THE 
C  EXACT  VALUE  OF  THE  QUANTITY  OF  AIR  ENTRAINED  IN  EACH  SEGMENT 

c********* *★*★★*★***** ***★★★★★★******★*************** ********** **** 

DO  42  KLM=1, 3 

0*  ***************************************************************** 

C  THIS  DO  LOOP  SOLVES  FOR  MOMENTUM  AND  ENERGY  CONSERVATION  IN  EACH 
C  OF  THE  EXISTING  SEGMENTS 

Q*  ***************************************************************** 

DO  100  IJ=IJl,J,2 

0*  ***************************************************************** 
C  ANALYSIS  OF  ODD  NUMBERED  SEGMENTS  IS  DONE  IN  THIS  SECTION 
0*  ***************************************************************** 
N=I 

JJ=IJ-1 

JJJ=IJ+1 

IF(VR1(JJ, I) . LT. 0. 1) VR1 ( JJ, I) =0 . 1 
IF(AVGD2(JJ,I) . LT. 1. E-4) AVGD2 ( JJ, I) =1. E-4 
IF ( I J.GT. 2) GO  TO  4000 

XRE1(JJ,I)=FXRE(VR1(JJ,I) ,AVGD2(JJ,I) ,C0NST1) 

GO  TO  4001 

4000  XRE1 (JJ, I) =FXRE (VRl (JJ , I— 1 ) , AVGD2 ( JJ , 1-1) ,CONSTl) 

4001  IF (XRE1 ( JJ, I) . LT . 1 2 . ) XRE1 ( J J , I ) =12. 

CD1 (JJ, I) =FCD(XRE1 (JJ, I) ) 

BCD (1,1) =CD1 ( JJ, I) 

IF ( I J .GT . 2) GO  TO  4007 

XC1(JJ,I)=FXC(AVGD2(JJ,I) ,CDl(JJ,I) ,C0NST2) 

GO  TO  4008 

4007  XC1 ( JJ , I) =FXC ( AVGD2 ( JJ , I— 1 ) ,CD1(JJ,I)  fC0NST2) 

4008  IF (XC1 ( JJ, I)  . LE.0.)XC1(JJ,I)=1. E-8 
C  WRITE ( 5 , * ) VRl ( JJ , I) 


IF(KLM.EQ.1)CALL  SWRL(P1,AVGD2(JJ, I) ,CD1(JJ,I) ,XTT(IJ) ,VL1 (JJ, I) , 
1RXL , THXL , VP) 

IF(IJ.GT.2)G0  TO  4002 

VR2 ( I J, I) =FVR(VR1 ( JJ , I) ,XT ( JJ) ,XT(IJ) ,XC1(JJ,I)) 

VR2 { IJ, I) =SQRT (VR2 ( IJ, I) **2.+VP**2.) 

GO  TO  4003 

VR2 (IJ, I) =FVR(VR1 (JJ, 1-1) ,XT ( JJ) ,XT(IJ) ,XC1(JJ,I)) 

VR2 (IJ, I) =SQRT (VR2 (IJ, I)  **2.+VP**2.) 

CONTINUE 

IF( I . EQ. 1) GO  TO  1200 
IF(IJ.EQ.2) GO  TO  1201 

TA ( IV, IJ-1 , I) = (FTA( IV, IJ-2 ,1-1) +TAC ( I V, 1 , 1-1) ) /2 . 

GO  TO  1205 

TA( IV, 1, I) =TAC ( IV, 1, 1-1) 

GO  TO  1205 
TA(IV, 1, 1) =T1 

IF(DS(JJ,N) .LE.l.E-4)GO  TO  61 
GO  TO  62 
DS (JJ,N)  =l.E-4 
CONTINUE 

IF ( I J. EQ. 2) GO  TO  1206 
TL=TBP 

IF ( I J.GT . 2) GO  TO  4004 

CALL  EVAP(TBP,P1,TA(IV, IJ-1, I) , HLF, SPHV, SHA) 

GO  TO  4005 

CALL  EVAP (TBP, Pl,TA ( IV, IJ-1, I) , HLF, SPHV, SHA) 

XKl ( JJ, I ) =FXK (XK00,XRE1 ( JJ, I) , SC) 

DS (IJ, I) =FD(DS (JJ ,N) ,XKl ( JJ , I) ,DT) 

IF(DS ( IJ, I) ) 31, 31,20 

DS  ( I  J,  I)  =1 .  E-8 

DS (I J, I) =SQRT (DS (IJ, I) ) 

AVGDS ( I J , I) =AVG (AVGDS (JJ,N) ,DS(IJ, I) ) 

IN=I 

IF(DS (IJ, I) -GT.DS(JJ,N) ) DS ( I J, I) =DS ( JJ , N) 

FV(IJ,I) =PI/6.*RHOL* (DS (JJ,N)  **3.-DS (IJ, I)  **3.)  *XNDS (I) 

FL  ( I  J,  I)  =FL(  JJ ,  N)  — FV  ( I  J,  I) 


C  ENERGY  BALANCE  STARTS 

Q*  *********************************  *  *****************  *  ************  * 

TAC(IV, IJ-1, 1) =TEMP (TA(IV, IJ-1, 1)  ,P1,P2) 

HLT=FV(IJ,I)*HLF 

TAM ( IV, IJ-1 , I) = (TA ( I V, IJ-1 , I) +TAC ( IV, IJ-1 , I) ) /2. 

HV=FV ( I J, I) *SPHV*TBP 

HA=XMAIR ( IV, IJ-1) *SHA*TAM( IV, IJ-1 , I) 

H1=HV+ HA-HLT 

H2=XMAIR(IV, IJ-1)*SHA+FV(IJ, I) *SPHV 
FTA( IV, IJ-1, I ) =Hl/H2 

q* ******  *******  ***********  ********  ************ *********************** 

C  ENERGY  BALANCE  IS  COMPLETE 

Q******************************************************************* 

TFV (IV, IJ-1, I) =TFV ( IV, IJ-2, 1-1) +FV ( IJ, I) 

IF(FL(IJ,I) .LE.0.)GO  TO  100 
VOLL ( I J , I ) =FL ( I J , I ) /RHOL 
D ( I J, I) = (6./PI/XND (I) *VOLL (I J , I))**. 333334 
IF(D ( I J, I ) >32,32,33 


32  D(IJ,I)=l.E-4 

33  AVGD2 (I J, I) =AVG (AVGD2 (JJ, 1-1) ,D(IJ,I)) 

35  N=I— 1 

JHA= (XT ( I J) -XT (JJ) ) /XC1 (JJ, I) 

IF(JHA. LT.0. 0001) JHA=. 0001 
VX(I  J,  I)  =VR2  (IJ,  I)  +XVAW 
VL1  (I  J,  I)  =VX  ( I  J,  I) 

IF ( I J. EQ. 2)  AVGVR2 ( IJ, I) =AVG (VR1 ( JJ , I) ,VR2(IJ,I)) 

IF(IJ.GT.2)  AVGVR2 (IJ, I) =AVG (VR1 ( JJ, 1-1)  ,VR2(IJ,I) ) 

VR1(IJ, I) =VR2(IJ, I) 

IF(AVGVR2(IJ, I) .LT. 0. 1) AVGVR2 ( I J, I) =0 . 1 
IF(AVGD2(IJ,I) .LT. l.E-4) AVGD2 ( IJ, I)  =l.E-4 
XRE2(IJ, I) =FXRE(AVGVR2 (IJ, I) ,AVGD2(IJ,I) ,C0NST1) 

IF(XRE2 ( I J , I )  .LE.12.) XRE2 ( IJ, I)  =12. 

CD2 ( I J, I) =FCD (XRE2 ( I J, I) ) 

CD1 ( I J , I ) =CD2 ( I J , I ) 

BCD (I , IJ) =CD1 (IJ, I) 

IF ( I J.GT. 2) GO  TO  4015 

AVGCD2 ( I J, I) =AVG (CD1 ( JJ , I) ,CD2(IJ,I) ) 

GO  TO  4016 

4015  AVGCD2 ( I J, I) =AVG (CD1 ( JJ , 1-1) ,CD2(IJ,I)) 

4016  XC2(IJ, I) =FXC(AVGD2(IJ, I) ,CD2(IJ,I) ,CONST2) 

IF (XC2 ( IJ, I) . LE.0.) XC2 (IJ, I) =l.E-8 

AVGXC  ( I J ,  I)  =AVG  (XC1  ( J  J ,  I)  ,  XC2  ( I J ,  I ) ) 

IF(IJ. GT.  I)  GO  TO  22 

C* ********************************************************** ******* 

C  EVEN  NUMBERED  SEGMENTS  ARE  ANALYSED  IN  THIS  PORTION  OF  THE  PROGRAM 
q*  ********************************** ******************************** 

IF(KLM.EQ. 1)  CALL  SWRL (P1,AVGD2 ( I J, I)  ,CDl(IJ,I)  ,XTT(JJJ) ,VL1(IJ,I)  , 
1RXL,THXL,VP) 

VR2 ( JJJ , I) =FVR (VR2 ( I J , 1-1) ,XT(IJ) ,XT(JJJ) ,XC2 ( I J, 1-1) ) 

VR2 ( J J J , I ) =SQRT ( VR2 ( J J J , I )  *  *  2 . +VP* *  2 . ) 

IF(DS(IJ,N) .LE. l.E-4) GO  TO  63 
GO  TO  64 

63  DS(IJ,N)=l.E-4 

64  CONTINUE 

TA ( IV, IJ, I) = (FTA ( IV, IJ-1,I-1)+TAC( IV, 1,1-1) )/2. 

CALL  EVAP (TBP, P1,TA( IV, I J, I)  , HLF,SPHV, SHA) 

XK2 ( I J, I) =FXK (XK00 ,XRE2 ( I J, 1-1) ,  SC) 

DS (JJJ, I) =FD (DS (I J,N) ,XK2 (IJ,I)  ,DT) 

IF(DS(JJJ, I))  34,34,21 

34  DS (JJJ , I) =l.E-8 

21  DS (JJJ , I) =SQRT (DS (JJJ , I) ) 

AVGDS (JJJ, I) =AVG (DS ( JJJ, I) ,AVGDS(IJ,N) ) 

IF  (DS  (JJJ ,  I)  .GT.DS  (I  J,  N) )  DS  (JJJ ,  I)=DS(IJ,N) 

FV (JJJ, I) =PI/6.*RHOL* (DS ( I J, N)  ** 3.-DS ( JJJ , I) ** 3 . ) *XNDS ( I) 

FL(JJJ, I) =FL(IJ,N) -FV( JJJ, I) 

Q*  ****************************************************************** 

C  ENERGY  BALANCE  : 

q** *********** ************************************* **************** 

TAC  ( IV,  I  J,  I)  =TEMP  (TA  ( IV,  I  J,  I)  ,P1,P2) 

HLT=FV (JJJ , I)  *HLF 

TAM( IV, IJ, I) = (TA( IV, IJ, I) +TAC( IV, IJ, I) ) /2. 

HV=FV (JJJ, I) *SPHV*TBP 

HA=XMAIR ( IV, IJ) *SHA*TAM( IV, IJ, I) 
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H1=HV+HA-HLT 

H2=XMAIR ( IV, I J) *SHA+FV ( JJJ , I) *SPHV 
FTA(IV, IJ, I) =H1/H2 

Q** ***********************************************************  ****** 

C  ENERGY  BALANCE  COMPLETE 

q************** ***********************************************  ****** 

TFV ( IV, I J, I) =TFV (IV, IJ-1, 1-1) +FV (JJJ, I) 

IF(FL(JJJ,I) .LE.0.)GO  TO  100 
VOLL ( J J J , I ) =FL ( J JJ , I ) /RHOL 
D (JJJ , I) = (6./PI/XND (I) *VOLL (JJJ , I) ) ** . 333334 
IF(D(JJJ,I)) 38,38,39 

38  D(JJJ, I) =l.E-4 

39  AVGD2 (JJJ , I) =AVG (D (JJJ , I) ,AVGD2(IJ,N) ) 

JHA1= (XT (JJJ) —XT ( IJ) )  /XC2 ( I J, I) 

IF(JHA1.LT.0. 000001) JHA1=. 000001 
AVGVR2 (JJJ ,1) =AVG (VR2 (I J, 1-1)  ,VR2(JJJ,I) ) 

IF ( AVGVR2 (JJJ , I) .LT. 0. 1) AVGVR2 (JJJ,I)=0.1 
IF(AVGD2 (JJJ, I) • LT. 1. E-4) A  VCD  2 (JJJ, I) =1.E~4 

XRE2 (JJJ , I) =FXRE(AVGVR2 (JJJ , I) ,AVGD2 ( JJJ, I) ,CONSTl) 

IF (XRE2 (JJJ, I)  .LE.12.) XRE2 (JJJ, I) =12. 

CD2 (JJJ, I) =FCD(XRE2 (JJJ, I) ) 

AVGCD2 (JJJ, I) =AVG(CD2 (IJ, 1-1) ,CD2(JJJ,I) ) 

VR1 (JJJ, I) =VR2 ( J J J , I ) 

VX(JJJ,I)  =VR2  (JJJ,  I)  +XVAW 
VL1 (JJJ , I) =VX (JJJ , I) 

22  CONTINUE 

100  CONTINUE 

XFV=0. 

DO  80  JKK=2, J 
XFV=XFV+  (FV (JKK, I) ) 

80  CONTINUE 

FVTOT=FVTOT+XFV 
SL=FL(1, I) 

IF ( I .EQ. 1) GO  TO  37 
DO  36  IDD=2, I 
IKK=I-1 

36  SL=SL+FL( IDD, IKK) 

37  TFL ( J, I) =SL-XFV 

TVOLL (J, I) =TFL(J, I) /RHOL 

C*** ********* ********************************************************* 

C  CONSERVATION  OF  MOMEOTUM:TO  SOLVE  FOR  THE  VELOCITY  OF  LEADING  DROPLET 

q************************** ******************************************** 


XM=0 . 

XN1=FL(1, I) *VL1(1, I) 

XFL=FL (1,1) 

XN2=0 . 

XN4=0. 

IF ( I . EQ. 1) GO  TO  50 
DO  25  ID=2, I 
I  K=  I D- 1 
IA=I+2-~ID 
IC=I+3-'ID 

XN1=XN1+FL( ID, I I) *VL1 (ID, II) 

XN2=XN2+FL (ID, I) *VL1 (ID, I) 

XN4=XN4+  (FV ( ID,  IK)  +XMAIR  ( IV,  IK) )  *XV AW 


LVU VMAT 


sera 


25  XFL=XFL+FL ( 1 , ID) 

50  CONTINUE 

XN3=FL(J, I) *VR2(J, I) 

XN=XNI -XN2-XN3 +XN4 
DN1=FL(J, I) 

DN2=0 . 

DO  26IM=2, J 
IKK=IM— 1 

26  DN2=DN2+FV (IM, IKK) +XMAIR(IV, IKK) 

DN=DN 1 +AG  *  DN2 

VAW=0 . 

VA=XN/DN 
VAW=AG*VA+VAW 
C  WRITE  ( 5 ,  * )  VAW 

IF(KLM.  LE.  3)  XXVAW=VAW 
VL1  (J,  I)  =VR2  ( J,  I)  +XXVAW 
RAJU=RAJ 

RAJU=RAJU*PI/180 . 

IF ( I .EQ.l) GO  TO  952 
GO  TO  953 

952  AVGVL=0 . 5* (VL1 (1,1) +VL1 (2,1)) 

GO  TO  954 

953  AVGVL=0. 5*  (VL1  (1 , 1-1) +VL1 (J, I) ) 

954  XT(J)= (XTP(I) +AVGVL*DT) 

XTT ( J) =XT ( J) *COS (RAJU) 

CALL  SWRL(P1,AVGD2(J, I) ,AVGCD2(J,I) ,XTT(J) ,AVGVL, RXL, THXL, VP) 

C****** ************************************************ **************** 

C  TO  DETERMINE  THE  GEOMETRY  OF  SPRAY  AND  EVALUATE  SPRAY  TIP  PENETRATION 

Q** *********** ************************************************** ******* 

ARXL (J) =RXL 
IF ( 1—2) 7,8,9 

7  XTP ( J) =ARXL ( J) 

SI (J)=THXL*180./PI 
GO  TO  1007 

C  IF ( IV.GT. 2) GO  TO  1007 

C  XXTP(J) =XTP(J) *COSD(SI (J) ) 

C  YXTP (J) =XTP (J) *SIND (SI (J) ) 

C  GO  TO  6 

8  SI GMA=PI-THXL 
VI J=COS (SIGMA) 

XTP ( J) = ( (XTP(J-l) ) **2.+ (ARXL(J) ) **2.-2.*XTP(J-l) *ARXL(J) *VIJ) ** .5 
GO  TO  10 

9  BETA1=- (XTP(J-2) )**2.+ (ARXL ( J— 1 ) ) **2 . + (XTP ( J — 1) ) **2 . 
BETA2=2.*ARXL(J-1) *XTP(J-1) 

BBETA=BETA1/BETA2 

IF(BBETA.EQ. 1.) BBETA=0. 9999999 
BETA=ACOS ( BBETA) 

BTA=BETA+THXL 
SV=PI-BTA 
AI=COS (SV) 

XTP(J) = ( (XTP(J-l) ) **2.+ (ARXL(J) ) **2.-2.*XTP(J-l) *ARXL(J) *AI) **0. 5 

10  XLAMDA1=- (ARXL (J) ) **2.+  (XTP(J-l) ) **2.+  (XTP (J)  )  **2. 
XLAMDA2=2.*XTP(J-1) *XTP(J) 

XXLAMDA=XLAMDA1/XLAMDA2 

I F ( XXLAMDA . EQ. 1 . ) XXLAMDA=0 . 99 9 99 9 9 


XLAMDA=ACOS (XXLAMDA) 

SI  (J)  =SI  (J-'l)  +XLAMDA*180./PI 
C  IF(IV.GT.2)G0  TO  1007 
C  XXTP(J) =XTP(J) *COSD(SI (J) ) 

C  YXTP (J) =XTP (J) *SIND (SI (J) ) 

C  GO  TO  6 

1007  XXTP( J) =XTP (J) *COSD ( (TH (IV)*180./PI)-~SI(J) ) 

YXTP ( J) =XTP (J) *SIND(- (TH (IV) *180./PI)+SI (J) ) 

6  WA=820. 

XTV=XXVAW*DT+XTV 
XLV =WA*  XTV*  DT 
THXV=ATAN (XLV/XTV) 

XTV=SQRT ( XTV* *  2 . +XLV** 2 ) 

ATHXV (J) =ATHXV ( I ) +THXV 
FTHXV(J) = ATHXV ( J) *180./PI 
XTPV (J) =XTV 

XXTPV ( J) =XTPV ( J) *COSD (FTHXV (J) + (TH ( IV) *180 ./PI) ) 

YXTPV ( J) =XTPV ( J) *SIND (FTHXV (J) + (TH (IV) *180. /PI) ) 
IF(KLM.GE.3)GO  TO  42 

*********************** 

C  TO  CALCULATE  ACTUAL  MASS  OF  AIR  ENTRAINED  IN  DEFORMED  CONE 

************************** ******************** ****************** 

CALL  CONEV (XTP ( J) ,TH(IV) ,XXTP(J) ,YXTP(J) ,XXTPV(J) ,YXTPV(J) , 
LXTP (J-'l)  ,DVOL) 

IF ( I . EQ. 1) GO  TO  1000 
VOA ( J ) =VOA ( I ) +DVOL 
GO  TO  1001 

1000  VOA ( J ) =VOA ( I ) +DVOL*2 ,/3 . 

1001  IF (IV.EQ.2) GO  TO  1004 

Q*  ****************************************************** *********** 

C  TO  CALCULATE  MASS  OF  AIR  IN  THE  ANNULAR  SECTION  OF  OUTER  CONES  BY 
C  INTERPOLATING  THE  VOLUME  OF  INNER  CONE  TO  BE  SUBTRACTED 

Q*  ************ **************************************************** 

CALL  CRAFT1 ( LK+1 , DUMMY2 , DUMMYl , XTP ( J ) ,FCOVOLI) 

GO  TO  1003 
1004  FCOVOLI=0.0 
1003  TCOVOL=VOA  ( J)  -'FCOVOLI 

I F ( TCOVOL . LE . 0 . ) TCOVOL= .001 
A VOA ( J ) =TCOVOL-TVOLL ( J , I ) 

XMAIR( IV, I)  =  (AVOA( J) -AVOA(I) )  *RHOA 
42  CONTINUE 

AREA(l) =AREA(2) 

AA(I) =AREA(2) 

XVAW=VAW 
XXVAW=VAW 
30  CONTINUE 

V(I)  =VAW 
X(I,  I)=0. 

J=I  +  1 

XI  ( I ,  I)  =X  ( I ,  I)  +XT  ( J) 

LLL=I 

XXFV(LLL) =XFV/XT ( J) 

YXT (LLL) =XT (J) -XT ( J-l) 

TIM(I) =TIME 
P1=P2 


V1=V2 

110  CONTINUE 

DO  1100  1=1, LK 
J=I+1 

C  WRITE (5 ,*)  XTP(J) 

C  WRITE (1,*) XTP(J) ,XXTP(J) ,YXTP(J) ,XXTPV(J) ,YXTPV(J) 

1100  CONTINUE 

DO  1008  J=2,LK+1 
DUMMY1 (J) =VOA(J) 

DUMMY2(J) =XTP(J) 

PEN (IV/J) =XTP (J) 

IF ( IV. EQ. 2) WRITE (5,*) PEN(IV,J) 

1008  CONTINUE 

DO  172  J=2,81,5 

C  WRITE (5 ,*) XXTP (J) ,YXTP(J) fXXTPV(J) ,YXTPV(J) 

172  CONTINUE 

111  CONTINUE 

C*********************************** ******************************* 

C  EVALUATION  OF  EQUIVALENCE  RATIO  WITHOUT  CONSIDERING  OVERLAP 

0t  ***************** ****** ****** *********************************** 

DO  4020  MN=2,5 
DO  4022  LN=1, LK 
DO  4021  NM=1,LN 

PHE (MN, NM) =TFV(MN, NM,  LN) / (XMAIR (MN, NM)  *0 . 06667) 

C  WRITE (5,*) PEN(MN,NM+1) , XMAIR (MN,NM) , TFV (MN, NM, LN) , PHE(MN,NM) 

4021  CONTINUE 

4022  CONTINUE 
4020  CONTINUE 

C***************************************************************** 

C  NCW  CONSIDERING  OVERLAP 

Q****************************** *********************************** 

DO  3000  J=1,LK+1 

DO  3001  L=2, 5 

DO  3002  IV=L+1, 5 

ADD3=0.0 

ADD5=0.0 

K=J 

3003  IF  (PEN  (L,  J)  -'PEN  ( IV,  K) )  3004,3005,3005 

3005  IF(K.GE.LK+1) GO  TO  3015 
K=K+1 

GO  TO  3003 

3004  M=K-1 
N=K+1 

3006  IF(PEN (L, J+l) -PEN ( IV, N) ) 3007,3008,3008 
3008  I F ( N . GE . LK+ 1 ) GO  TO  3007 

N=N+1 

GO  TO  3006 

3007  IJ=N-1 

Dl= (PEN (L, J) -PEN ( IV,M+1) ) / (PEN ( IV,M+1) -PEN ( IV, M)  )  *TFV ( IV,M, LK) 
ADD1=ABS (Dl) 

D2= (PEN (L, J+l) -PEN ( IV, IJ) )/ (PEN (IV, I J+l) -PEN (IV, IJ))*TFV(IV, IJ 
ADD2=ABS(D2) 

IF (N.GE. LK+1) ADD2=0 . 0 
IF ( I J. GE.M+2) GO  TO  3009 

IF( I J. EQ. 60 .AND. M, EQ. 59) ADD5=TFV (IV, IJ,LK) 
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GO  TO  3010 

3009  DO  3011  LM=M+2,IJ 
ADD3=TFV(IV,LM~1,LK) 

3011  CONTINUE 

3010  ADD4=ADD1+ADD2+ADD3+ADD5 

XNV= (IJ+M) /2.*0. 31+ (TH{IV) *180. /PI) 

XDE=J*0. 31+ (TH(L) *180. /PI) 

RATIO= (TAND (XDE) ) **2./ (TAND (XNV) )**2. 

ADD=RATIO*ADD4 

TFV ( L, J , LK) =TFV ( L , J , LK) +ADD 

TFV (IV,M, LK) =TFV ( I V,  M,  LK) -RATIO* ADD1 

TFV ( IV, I J, LK) =TFV ( IV, I J, LK) -RATIO* ADD2 

IF ( I J. LT.M+2) GO  TO  3002 

DO  3013  KL=M+2, IJ 

TFV (IV,KL-1,LK) =TFV ( IV, KL-1, LK) -RATIO* (ADD 3/ ( I J-M-l) ) 

3013  CONTINUE 
3002  CONTINUE 

3015  IF(XMAIR (L, J) .EQ.0.0)GO  TO  3001 

PHE (L, J)  =TFV (L, J,LK) / (XMAIR (L,J) *0.06667) 

DELAY (L,J,LK) =A* (P1**AN) * (PHID**C) * (EXP (DD/FTA (L, J, LK) ) ) 

WRITE (5,*) ' PHE’ , PHE (L, J) , ’DELAY' , DELAY (L,J,LK) , 'J' ,J 

IF(DELAY(L, J,LK) .GT. (J*DT*1000.) )GO  TO  3001 

I=LK 

GO  TO  5001 
3001  CONTINUE 
3000  CONTINUE 

5000  CONTINUE 

5001  NO=0 

0** ****************************************************  *********** 

C  TO  CALCULATE  THE  PRESSURE  RISE 

0*  **************************************************************** 

SMRT=0 
SMRT1=0 
DO  5002  K=2, 5 
DO  5003  L=l, I 

IF ( (PHE (K, L) .GT.0.3) .AND. (PHE(K,L) .LT.3.) ) GO  TO  5004 
SMRT1=SMRT1+ (TFV { K, L, I ) +XMAIR (K,L) ) *2. 83*FTA(K, L, I) 

SAIR=SAIR+ XMAIR (K,  L) 

GO  TO  5003 

5004  CALL  AFT(PHE(K,L)  ,FTA(K,L, I) ,Pl,AT,ARR) 

SMRT=SMRT+ ( TFV ( K, L , I ) +XMA IR ( K , L) ) * ARR* AT 
SAIR=SAIR+XMAIR (K, L) 

NO=NO+l 

5003  CONTINUE 

5002  CONTINUE 

C  SINCE  THERE  ARE  4  NOZZLES 
SMRT 1 = 4  *  SMRT 1 

SMRT=4*SMRT*2. 2*778*12*2. 54**3/14700 
C  IF  SPIKE  IS  ASSUMED  TO  BE  3  DEGREES  CRANK  ANGLE 
C  THE  VOLUME  AND  TEMPERATURE  DUE  TO  COMPRESSION  ALONE  WOULD  BE: 
ANG=ANG-3* PI/180. 

CALL  PV(ANG,V1, PI, V2, P2) 

TAIR=TEMP(TA(2,1,I)  ,Pl,P2) 

Ell=Vl*RHOA-4*SAIR 

E1=E11*2.83*TAIR 


Pl= (SMRT1+SMRT+E1) /V2 
C  WRITE ( 1 ,  * ) ANG* 180/PI , PI 

V1=V2 

C  MASS  OF  FUEL  COMING  OUT  IN  THE  REMAINING  TIME  FROM  ALL  4  NOZZLES 
RM=4* (PI/4 .*DNN*DNN*VLi (I, I) *RHOL*3/(6*RPM) ) 

C  ASSUMING  THIS  FUEL  TO  BURN  AT  PHI=1  IN  NEXT  1  DEGREE  OF  ROTATION 
TT=TAIR 

CALL  AFT (1. ,TAIR, Pl/AT, ARR) 

SMRT=SMRT+ ( RM+RM) *ARR*AT*2. 2*778*12*2. 54**3/14700 
C  BECAUSE  PHI=1  MASS  OF  AIR  WILL  ALSO  BE  RM 
ANG=ANG— PI/180. 

CALL  PV (ANG, VI, PI, V2, P2) 

TAIR1=TE*1P  (TT,  PI,  P2) 

El= (Ell-RM) *2. 83*TAIR1 
P2= (SMRT1+SMRT+E1)  /V2 
C  WRITE (1,*)ANG*180/PI,P2 

c* ****★*★★★*★**★*★★**★★★★**★★★★*★★**★★★*★*★***★*★*★*★★★★**★*★*★**★**** 

C  THIS  IS  A  PARTICULAR  SOLUTION  AND  WILL  REQUIRE  CHANGES  FOR  EACH  TEST 

C  RUN.  IT  IS  POSSIBLE  TO  SUBSTITUTE  IT  WITH  A  GENERAL  SOLUTION 

£***★* a*************************************************************** 

5010  P1=P2 
V1=V2 

ANG=ANG-PI/180 . 

CALL  PV (ANG, VI, P1/V2, P2) 

C  WRITE (1,*)ANG*180/PI,P2 
IF  (P2. LT. 5) GO  TO  5011 
GO  TO  5010 

5011  CLOSE  (UNIT=1) 

STOP 

END 

Q*** *********************************************************** ******** 

SUBROUTINE  SWRL (PT,D,CD,XT,AVGVL,RXL,THXL,VP) 

DOUBLE  PRECISION  RXL, THXL,WA1,WA2,WL,XL 

PI-3.1416 

GC-981. 

TAO=305. 

PTO=l. 


R=2.8326 

RHOAO=PTO/(R*TAO) 

THINJ=11 

IM12=60 

DELAY=8. 

RPM=900. 

DT=THINJ/ ( IM12*6.*RPM) 

RHOA= (PT/PTO) **  (1./1.256) *RHOAO 

RHOL=. 74871 

CONST2=RHOL/RHOA 

WA=820. 

I F ( XT . LT . 0 . 01 ) XT=0 .01 
WA1=4  .*CONST2*GC*D/  ( 3  .*CD*XT*DT) 
WA2=SQRT ( (WA1+2 . * WA)  ** 2 .  -4  . *WA*WA) 
WL= (WA1+2. * WA-WA2 ) / 2 . 

IF ( AVGVL . EQ. 0 .) AVGVL=0 . 00001 

VP=AVGVL*DT*WL 

XL=VP*DT 
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RXL=SQRT ( (AVGVL*DT) **2.+XL*XL) 

THXL=ATAN (XL/ (AVGVL*DT) ) 

RETURN 

END 

q*  **★******★★★***★*****★**★★*★★★★★★******★****★★**★★★★★★*★*★****★★★*★★* 

SUBROUT INE  CONEV (X2 , ANGLE , XL, YL , XV, YV, XI , DV) 

COMMON/ A/ AREA ( 2 ) 

PI=3. 1415 

AMINOR=2*X2*TAN (ANGLE) 

AMAJOR= ( (XL-XV) **2.+ (YL-YV) **2.) **0.5 
AREA ( 2) =PI*AMINOR* AMAJOR 
BASE= (AREA (2) +AREA (1) )  /2. 

HT=X2-X1 

DV=BASE*HT 

RETURN 

END 

Q*************************************************** ******************** 

SUBROUTINE  CRAFT1 (MM,A,F,X,VRS) 

DIMENSION  F(61) ,D(61) ,B{61) ,DD(62) ,XF(0:61,61) 

DOUBLE  PRECISION  A(61),X 
DO  20  J=1 ,MM 
D (J) =A(J) -X 
IF  (J.EQ.l)GO  TO  40 

IF(ABS(D(J) ) .LT.ABS (D (J-l) ) -AND.D(J) .LT.0.0)GO  TO  40 
GO  TO  20 
40  K=J 

20  CONTINUE 

IF (K.EQ. 1) GO  TO  60 
DO  50  L=1,MM+1-K 
B(L)=A(K-1+L) 

XF(0,L) =F (K— 1+L) 

50  CONTINUE 

DO  70  M=MM+2-K,MM 
DO  70  N=1,K-1 
B(M) =A ( K-N ) 

XF ( 0 , M) =F(K-N) 

70  CONTINUE 

GO  TO  90 

60  DO  80  N=1,MM 

B(N)=A(N) 

XF(0/N)  =F(N) 

80  CONTINUE 

90  DO  100  JJ=1,MM 

DD ( JJ) =B ( JJ) -X 
100  CONTINUE 

DO  120  L=l,l 
DO  110  N=L+1, 2 

XF (L, N) = (XF (L-l, L) *DD (N) -XF (L-l, N) *DD(L) ) /(B(N) -B(L) ) 

110  CONTINUE 

120  CONTINUE 

VRS=XF(1,2) 

RETURN 

END 

Q*  ************************************************************************ 

SUBROUTINE  EVAP(TL,P,T,HLF,SPHV/SHA) 


-V 


C  TC  IS  CRITICAL  TEMP 

C  PC  IS  CRITICAL  PRESSURE 

C  XMW  IS  MOL.  WT.OF  FUEL 

FC=18. 

TC=660. 

R=1.987 

XMW=226. 

TM=  (T-TL) /ALOG (T/TL) 

HLF=(2.303*R*TL*TC*ALOG10(PC) )/( (TC-TL) *XMW) 
X=-9.G8+0.27045*TM-1.0133E-8*TM*m 
SPHV=(X+(5.03*P/PC*  (TC/™)  **3.) ) /XMW 
SHA=0. 24332+0. 00004593*™ 

RETURN 

END 

C******************************************************************* 

SUBROUTINE  PV (THETA, VOL1 , PR1, VOL2, PR2) 

C  THIS  SUBROUTINE  CALCULATES  VOLUME  OF  CYLINDER  FROM 
C  ENGINE  GEOMETRY. 

PI =3. 142 

q****************************************************************+ 

C  CL=CL£ARANCE,D=DIA  OF  CYLINDER, CRR=CRANK  RADI US , CR=LENGHT 
C  OF  CONNECTING  ROD,DRP=DIFF  BETWEEN  CENTERS  OF  CYLINDER 
C  AND  PISTON.  ALL  DIMENSIONS  ARE  IN  CM 

***★*★**★★****★★★***★**★★★*****★**★**★*★  ★★****★**★**•★**** ** 

CL=0 . 2841 
D=8.026 
CR=8 .89 

C  CRR  IS  APPROXMATE  VALUE 

CRR=2.29 
ANC=1. 39 

C  DRP  IS  APPROXIMATE  VALUE 

DRP-0.006 

Q******** ************ ************************* ********** *********** 

C  PR  IS  PRESSURE  VOL  IS  VOLUME  OF  CYLINDER 
£********************★********************************************* 

VOLll=SQRT (CR**2- (DRP+CRR*SIN (THETA) ) **2) 

VOL2=PI/4.* (D**2) * (CL+CR+CRR* (l.-COS (THETA) ) "VOL11) 

PR2=PR1* (VOL 1/VOL 2) **ANC 
C  WRITE ( 5 , * ) VOL2 , PR2 

RETURN 
END 

FUNCTION  TEMP (Tl, Pi, P2) 

TEMP=T1* (P2/P1) **0.286 
C  WRITE (5,*)  TEMP 

RETURN 
END 

C********************************************************************** 

SUBROUTINE  AFT(FF,TT,PP,AT,ARR) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/BLOCK/AN, AM, AL,AK,F,T,P,KL0,IERR,XPA 
COMMON/TAB/DXT (12) ,DXP(12) ,DXF(12) 

COMMON/ PROP/ AVM , R, H,U, DRT,DHT,DUT,DRP,DHP,DUP,DRF,DHF,DUF 
F=FF 


AL=0. 

AK=0 . 

P=14.7*PP 
TM=TT* 1.8 
TH1=537.0/1.8 

HRF= (3 . 55E-5*TT**2-0 . 076997*TT~69. 95) *1000. *453 . 59/252. 

HAIR*  (1.286E-5*TM**2+0.274675*'IM+3. 33)  *29.8 

XAIR=24 . 5*4 . 607/F 

HMR=226.+XAIR*29.8 

HR= (HRF+XAIR*29 . 8) /HMR 

KL0=1 

T=6000. 

IND=0 

10  IND=IND+1 

CALL  PER(1,RR) 

ARR=RR 

C  WRITE (1,*) 'RR' ,RR,'NO.  OF  ITERATIONS IND 

IF (IERR.NE.0)  GO  TO  20 
889  FORMAT ( IX , ' H  HR  DELT  DHT' ,4F15. 5) 

C888  FORMAT ( IX, 'H  HR  DHT',3F10.2) 

DELT= (H-HR) /DHT 

T=T~DELT 

AT=T 

C  WRITE ( 1 , 889) H, HR, DELT , DHT 

IF (DABS (DELT)  .LE.l.)  GO  TO  30 
IF(IND.LT. 100)  GO  TO  10 
C  WRITE (1,40) 

40  FORMAT (10X,' ITERATION  DID  NOT  CONVERGE  IN  25  ATTEMPTS’) 

GO  TO  5 

30  WRITE (5, 60) AT 

60  FORMAT (10X, ’FLAME  TEMPERATURE  (DEG.  R)  =’,F8.1) 

GO  TO  5 

20  WRITE(5,70) IERR, AN, AM, AL, AK,F,T,P,KL0 

70  FORMAT (10X, 'SUBROUTINE  PER  FAILED.  ERROR  CODE  IS',I2/10X 

C,  'CONDITIONS  OF  FAILURE  (AN,AM,AL, AK,F,T, P, KL0) ARE' ,7F10. 2 

D,  15) 

5  WRITE (5,*)  'ARR' ,ARR 

RETURN 
END 

SUBROUTINE  CQMD(JDR, RR) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/BLOCK/ AN ,AM,AL,AK,PHI,T,P,KL0, IERR , XI , X2 , X3 , X4 , X5 , 
E  X6,X7,X8,X9,X10,X11,X12,X13 

COMMON/TAB/DX1T ,  DX2T ,  DX3T ,  DX4T ,  DX5T ,  DX6T ,  DX7T ,  DX8T ,  DX9T , 
FDX10T,DX11T,DX12T,DX1P,DX2P,DX3P,DX4P,DX5P,DX6P,DX7P,DX8P, 
GDX9P, DX10P , DX1 IP , DX1 2P , DX1F , DX2F , DX3F , DX4F, DX5F , DX6F, DX7F , 
HDX8F,DX9F,DX10F  ,.DX11F,DX12F 
C  COMMON/TAB/DXT ( 12) ,DXP(12) ,DXF(12) 

DIMENSION  A(4,4) ,B(4) ,C(4,3) 

DATA  JF , PREC/ 0,1. 0E-0  3/ 

IF( (T-1080.) * (7200. -T) ) 102,105,105 
102  IERR=L 

GO  TO  710 

105  R0= (AN+0. 25*AM-0. 5*AL) /PHI 

R=R0+0.5*AL 
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R1=R0*3. 7274+0. 5*AK 
R2=R0*0. 0444 
GG=R+R1+R2 

C  WRITE (6 ,*) '  GG= ’ ,GG 

999  FORMAT (IX, 'PHI  R  R0  Rl  R2',5F10.5) 

8  FORMAT (3F8. 4) 

IF ( R. GT. .  5* AN)  GO  TO  110 
IERR=2 
GO  TO  710 
110  Dl=AfVAN 

D2=2.0*R/AN 
D3=2. 0*R1/AN 
D4=R2/AN 

115  SQP=DSQRT (P/14.696) 

TA=0.005*T/9.0 
ALTA=DLOG (TA) 

TAIN=1. 0/TA 
TASQ=TA*TA 

998  F0RMAT(1X,'I  AM  O.K. ') 

Cl= (10. 0**(0.432168*ALTA-11.2464*TAIN+2. 67269-0. 745744E-01*TA 
A+0. 242484E-'02*TASQ) )  /SQP 

C2= (10.0** (0. 31 080 5*ALTA-12. 954 0*TAIN+3. 21779-0. 73 833 6E~01*TA 
B+0.344645E-02*TASQ)  )/SQP 

C3= (10. 0**(0.389716*ALTA-24.5828*TAIN+3. 14505-0. 963730E~01*TA 
C+0.585643E-02*TASQ) )/SQP 

C5=10. 0** (-0 . 14 1784* ALTA-2. 13308*TAIN+0 . 853461+0 . 355015E-01*TA 
D-0. 310227E-02*TASQ) 

C7=10.0** (0. 150879E-01*ALTA-4 . 70959*TAIN+0 .646096+0 . 272805E-02*T 
EA-0. 154444E-02*TASQ) 

C9= (10.0** (-0. 752364*ALTA+12. 4210*TAIN-2. 60286+0 . 259556*TA 
F-0.162687E-01*TASQ) ) *SQP 

C10= (10 . 0** (-0 . 415302E-02*ALTA+14 . 8627*TAIN-4 . 75746+0 . 124699*TA 
G-0.900227E-02*TASQ) ) *SQP 
C 

C  SECTION  DECIDES  WHETHER  TO  MAKE  A  NEW  ESTIMATE  OF  X4,X6,X8,X11 

C 

IF(KL0-1) 305,205,410 
205  IF(JF.EQ.0) GO  TO  305 

IF(DABS (PHI-PHIPR) . GT. 1.0E-6)GO  TO  305 
IF (DABS (T/TPR-1. 0) .GT.0.02)  GOTO  305 
IF(DABS (P/PPR-1 . 0) .GT. 0.05)GO  TO  305 
X4=X4PR 
X6=X6PR 

IF(X6.LE.0.0)X6=0. 0000001 
X8=X8PR 
X11=X11PR 
GO  TO  410 
C 

C  SECTION  300  CAN  MAKE  AN  INITIAL  ESTIMATE  OF  X4,X6,X8,_X11 

C 

305  IF(PHI.GE.1.0) GO  TO  310 

PAR=1. 0/ (R1+R2+R+0. 25* AM) 

GO  TO  315 

310  PAR=1.0/  (R1+R2+A.  ,+u.5*AM) 

315  FUN1=2.0*AN*C10 


FUN2=0.5*AM*C9 
FUN3=2. 0/PAR 
FUN4=2.0*R 
OX=1.0 

SQOX=DSQRT (OX) 

FOX= (FUNl*SQOX+AN) / (I . 0+C10*SQOX) +FUN2*SQOX/ ( 1 . 0+C9*SQOX) 
A+FUN3*OX~FUN4 

IF(FOX) 325,330,335 
OX=0.1*OX 

IF (OX.GE. 1. 0E-30) GO  TO  320 

IERR=3 

GO  TO  710 

IND=1 

FOX= ( FUN I *  S  QOX+AN) / ( 1 . 0+C10*SQOX) +FUN2*SQOX/ 

1  (1. 0+C9*SQOX) +FUN3*OX-~FUN4 

DOX=0 . 25*FUN1/(SQ0X* ( 1 . 0+C10*SQOX) **2) +0 . 5*FUN2/ ( SQOX* 

2  (1.0+C9*SQOX) **2) +FUN3 
RAT=FOX/DOX 
OX=OX~RAT 
SQOX=DSQRT (OX) 

IF (DABS (RAT/OX) .LE. 1. 0E-2) GO  TO  330 
IND=IND+1 

IF( IND. LE. 100) GO  TO  327 
X4=0.5*AM*PAR/(1.0+C9*SQOX) 

X6=AN*PAR/ ( 1 . 0+C 10*SQOX) 

IF  (X6.LE.0.0)X6=0. 0000001 
X8=OX 

Xll=Rl*PAR 


SECTION  400  CALCULATES  THE  ELEMENTS  OF  THE  MATRIX (LINEARISED  EQ) 
IND=1 
NCALL=0 

IF(X4 . LE. 0. 0)  X4=0. 0000001 
SQX4=DSQRT(X4) 

IF  (X8.LE.0.0)  X8=0. 0000001 
SQX8=DSQRT (X8) 

SQX11=DSQRT (Xll) 

Xl=Cl*SQX4 

X2=C2*SQX8 

X3=C3*SQX11 

X5=C5*SQX4*SQX8 

X7=C7*SQX11*SQX8 

X9=C9*X4*SQX8 

X10=C10*X6*SQX8 

T14=0. 5*C1/SQX4 

T28  =0 . 5*C2/ SQX8 

T311=0. 5*C3/SQX11 

T54=0. 5*C5*SQX8/SQX4 

T58=0 . 5*C5*SQX4/SQX8 

T78=0. 5*C7*SQX11/SQX8 

T711=0 . 5*C7*SQX8/ SQX11 

T94=C9*SQX8 

T98=0. 5*C9*X4/SQX8 

T106=C10*SQX8 

T108=0. 5*C10*X6/SQX8 
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A( 1, 1)  =T14  +  2.+T54+2.*T94 
A(l, 2) =-Dl* (1.+T106) 

A( 1, 3) = (T58+2.0*T98) -Dl*T108 
A(l,4) =0.0 
A(2, 1) =T54+T94 

A(2f  2)  =  (1. +2.*T106) -D2* (1.+T106) 

A(2,3) = (T28+T58+T78+2.+T98+2.*T108) -D2*T108 
A(2, 4) =T711 
A(3,l) =0. 

A  (3, 2) =-D3* (1.+T106) 

A(3/3) =T78~D3*T108 
A(3,4) =T311+T711+2. 

A (4,1)  =T14+1.+T54+T94 
A(4, 2) =1.+T106+D4* (1.+T106) 

A(4,3) =T28+T58+T78+1.0+T98+T108+D4*T108 
A(4,4)=T311+T711+1. 

IF(NCALL.EQ.l) GO  TO  810 
B(l)=- (X1+2.0*X4+X5+2.0*X9) +D1* (X6+X10) 

B (2) =- (X2+X5+X6+X7+2.*X8+X9+2.*X10) +D2* (X6+X10) 

B(3) =- (X3+X7+2.*X11) +D3* (X6+X10) 

B(4)  (X1+X2+X3+X4+X5+X6+X7+X8+X9+X10+X11+D4* (X6+X10) ) +1. 

SECTION  500:MATRIX  IS  SOLVED  BY  GAUSSIAN  ELIMINATION 

DO  505  K=l, 3 
KPl=K+l 

BIG=DABS (A(K, K)  ) 

IF(BIG.GE. 1. 0E-05)  GO  TO  520 
IBIG=K 

DO  510  I=KP1,4 

IF(DABS (A( I , K) ) .LE.BIG)  GO  TO  510 
BIG=DABS (A( I,K) ) 

IBIG=I 
510  CONTINUE 

IF (BIG.GT. 0. )  GO  TO  512 

IERR-4 

GO  TO  710 

512  IF ( IBIG. EQ.K)  GO  TO  520 

DO  515  J=K, 4 
TERM=A(K, J) 

A(K, J) =A( IBIG, J) 

A ( IBIG, J) =TERM 
515  CONTINUE 

TERM=B(K) 

B (K) =B ( IBIG) 

B(IBIG) =TERM 
520  DO  525  I=KPl, 4 

TERM=A  ( I ,  K)  /A(K,  K) 

DO  530  J=KP1,4 
A(I,J)=A(I,J)-A(K,J) *TERM 
530  COOTINUE 

B ( I) =B ( I ) -B (K) *TERM 
525  CONTINUE 

505  COETTINUE 

IF (DABS (A(4,4) ) .GT.0.)  GO  TO  550 
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IERR=4 
GO  TO  710 

550  B (4) =B (4) /A(4,4) 

C  WRITE (6,698) B(3) ,X3f X7,X11,X6,X10,D3 

698  FORMAT ( IX,’ B  X3  X7  Xll  X6  X10  D3',7F10.3) 

B(3)=(B(3)~A(3,4) *B(4) )/A(3,3) 
B(2)=(B(2)-A(2,3)*B(3)-A(2,4)*B(4))/A(2,2) 

B(l)  =  (B(1)~A(1,2)  *B(2)-A{  1, 3)  *B(3)-~A(1, 4)  *B(4)  )/A(l,  1) 

C 

C  CHECKS  PRECISION 

C 

602  NCK=0 

X4=X4+B(1) 

IF (DABS (B( 1) /X4  ) .GT.PRBC)  NCK=1 
X6=X6+B(2) 

IF  (X6.LE.0.0)X6=0. 0000001 
IF(DABS(B(2)/X6  ) .GT.PRBC)  MCK=1 
C  WRITE(1,699)X8,B(3) ,T311,T711 

699  FORMAT (IX, 'X8  B{3)  T311  T711* ,4F15.5) 

X8=X8+B(3) 

IF(DABS(B(3)/X8  ) .GT.PRBC)  NCK=1 
X11=X11+B(4) 

IF (DABS (B(4) /Xll) .GT.PRBC)  NCK=1 
XABS4=DABS(B(1)/X4) 

XABS6=DABS(B(2)/X6) 

XABS8=DABS(B(3)/X8) 

XABS 11=DABS ( B ( 4 ) /Xll ) 

C  WRITE (5,*) XABS4,XABS6,XABS8,XABS11 

C  WRITE (5,*) 'Bl' ,B(1) ,'B2' ,B(2) ,'B3' ,B(3) , 'B4 ' ,B(4) , ' NCK' ,NCK 

995  FORMAT ( IX,' X4  X6  X8  X11',4F15.5) 

C  WRITE ( 1 , 995 ) X4 , X6 , X8 , XI 1 

IF  (X8.LE.0.0)  X8=0. 0000001 
IF  (X4.LE.0.0)  X4=0. 0000001 

IF( (X4. GT.0.)  .AND. (X8.GT.0.) .AND. (Xll. GT.0.) )  GOTO  620 
IERR=5 
GO  TO  710 
C 

C  SECTION  700  CALCULATES  THE  REMAINING  MOLE  FRACTIONS 

C 

991  FORMAT ( IX, 'NCK=' ,12) 

620  IF(NCK.EQ.0)  GO  TO  625 

IF(IND.LT. 100)  GO  TO  622 
IERR=6 
GO  TO  710 
622  IND=IND+1 

GO  TO  455 

625  IF(X6.GE.0.)  GO  TO  702 

IERR=5 

997  FORMAT (IX, ' I  AM  ALL  RIGHT') 

GO  TO  710 
702  IERR=0 

JF=1 

PHIPR=PHI 

TPR=T 

PPR=P 


87 


* n V  W V 


I1/*1  T.1!*!  I V»  ' 


X4PR=X4 
X6PR=X6 
X8PR=X8 
X11PR=X11 
SQX4=DSQRT (X4) 

SQX8=DSQRT(X8) 

SQX11=DSQRT (Xll) 

X1=C1*SQX4 

X2=C2*SQX8 

X3=C3*SQX11 

X5=C5*SQX4*SQX8 

X7=C7*SQX11*SQX8 

X9  =C9  *X4  *  SQX8 

X10=C10*X6*SQX8 

X13= (X6+X10) /AN 

X12=R2*X13 

C  WRITE (I,*) '  X1’,X1,'X2',X2,'X3',X3,'X4’,X4,’X5',X5,'X6',X6, 

C  1  'X7' , X7 , ' X8 ' ,X8, 'X9' ,X9, 'X10' ,X10, 'Xll' ,X11, 'X12' ,X12, 'X13' ,X13 
IF(JDR-l) 725,800,800 
725  RETURN 

710  JF=0 

RETURN 

C 

C  SECTION  800  CALCULATES  THE  ELEMENTS  OF  THE  MATRIX  OF  THE 

C  PARTIAL  DIFFERNTIAL  EQUATIONS 

C 

800  NCALL=1 

GO  TO  460 

810  ZZ=ALOG(10.)  *0.005/9. 

DC1T=X1*ZZ* (0 . 432168*TAIN+11 . 2464/TAS(>-0 . 745744E-1+. 484968E-2*TA) 

DC2T=X2*ZZ* ( . 310805*TAIN+12. 9540/TASQ-. 738336E-1+.689290E~2*TA) 

DC3T=X3*ZZ* ( . 389716*TAIN+24 . 5828/TASQ-. 963730E-1+1 . 171286E-2*TA) 

DC5T=X5*ZZ* (-.141784*TAIN+2.13308/TASQ+.355015E-~1-.620454E— 2*TA) 

DC7T=X7*ZZ* ( . 150879E-l*TAIN+4. 70959/TASQf. 272805E-2-. 308888E-2*TA) 

DC9T=X9*ZZ* (-'.752364*TAIN-12.4210/TASQ+. 259556". 325374E— 1*TA) 

DC10T=X10*ZZ* (-. 415302E-2*TAIN~14 . 8627/TASQf . 124699-. 01800454*TA) 

PP2= . 5/P 

DC1P=-X1*PP2 

DC2P=-X2*PP2 

DC3P=-X3*PP2 

DC9P=X9*PP2 

DC10P=X10*PP2 

D5AN=-R0*X13/PHI 

DD4F=0.0444*D5AN 

C ( 1 , 1) =- (DC1T+DC5T+2 . *DC9T-D1*DC10T) 

C(2,l) =- (DC2T+DC5T+DC7T+DC9T+ (2.-D2) *DC10T) 

C  ( 3 , 1) =- (DC3T+DC7T-D3*DC10T) 

C ( 4 , 1 ) =- (DC1T+DC2T+DC3T+DC5T+DC7T+DC9T+ ( 1. +D4) *DC10T) 

C(l,2)=— (DCIP+2.*DC9P-D1*DC10P) 

C (2, 2) — (DC2P+DC9P+ ( 2. ~D2) *DC10P) 

C(3, 2) =- (DC3P-D3*DC10P) 

C(4,2) =- (DC1P+DC2P+DC3P+DC9P+ (1.+D4) *DC10P) 

C(l,3)=0. 

C(2,3) =2.*D5AN 
C(3,3) =7.4548*D5AN 
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C(4,3)=-DD4F 
DO  905  K=l,3 
KP1=K+1 

AMAX=DABS (A(K,K) ) 

MAX=K 

DO  910  I=KP1,4 

IF(DABS(A(I,K) ) .LE.AMAX)  GO  TO  910 
AMAX=DABS ( A  ( I ,  K) ) 

MAX=I 

910  CONTINUE 

IF(AMAX. GT. 0.0)  GO  TO  912 

IERR=7 

GO  TO  710 

912  IF(MAX.EQ.K)  GO  TO  950 

DO  915  J=K,4 
TEEM=A(K,J) 

A(K, J) =A(MAX,J) 

A(MAX,J)=TERM 
915  CONTINUE 

DO  920  J=l,3 
TEPM=C(K,J) 

C(K,J)=C(MAX,J) 

C(MAX,J)=TERM 
920  CONTINUE 

950  DO  925  I=KP1,4 

TERM=A(I,K) /A(K,K) 

DO  930  J=KPl/4 
A( I, J) =A( I, J) — A(K, J)  *TERM 
930  CONTINUE 

DO  935  J=l, 3 

C ( I , J) =C( I , J) ~C(K, J) *TERM 
935  CONTINUE 

925  CONTINUE 

905  CONTINUE 

IF(DABS (A(4,4) ) .GT.0.0)GO  TO  938 

IERR=7 

GO  TO  710 

938  DO  940  J=l,3 

C(4,  J)  =C(4,  J)/A(4 ,4) 

C(3,J)=(C(3,J)-A(3,4)*C(4,J))/A(3,3) 

C(2,J)  =  (C(2,J)-A(2,3)  *C(3,J)-A(2,4)  *C(4,  J)  )/A(2,2) 
C(1,J)=(C(1,J)-A(1,2)*C(2/J)-A(1,3)*C(3,J)-A(1,4)*C(4,J))/A(1,1) 
940  CONTINUE 

DX4T=C(1,1) 

DX6T=C(2,1) 

DX8T=C(3,1) 

DX11T=C(4,1) 

DX4P=C(1,2) 

DX6P=C(2,2) 

DX8P=C(3,2) 

DX11P=C(4,2) 

DX4F=C(1,3) 

DX6F=C(2,3) 

DX8F=C(3,3) 

DX11F=C(4,3) 
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DX1T=T 14  *DX4T+DC1T 

DX2T=T28*DX8T+DC2T 

DX3T-T311*DX11T+DC3T 

DX5T=T54*DX4T+T58*DX8T+DC5T 

DX7T=T78*DX8T+T711*DX11T+DC7T 

DX9T=T94*DX4T+T98*DX8T+DC9T 

DX10T=T106*DX6T+T108*DX8T+DC10T 

DX12T=D4* (DX6T+DX10T) 

DX1P=T14*DX4P+DC1P 
DX2P=T28*DX8P+DC2P 
DX3P=T311*DX11P+DC3P 
DX5P=T54  *DX4P+T58  *DX8  P 
DX7P=T78*DX8P+T711*DX11P 
DX9P=T94*DX4P+T98*DX8P+DC9P 
DX10P=T106*DX6P+T108*DX8P+DC10P 
DX12P=D4* (DX6P+DX10P) 

DX1F=T14*DX4F 

DX2F=T28*DX8F 

DX3F=T311*DX11F 

DX5F=T54*DX4F+T58*DX8F 

DX7F=T78*DX8F+T711*DX11F 

DX9F=T94*DX4F+T98*DX8F 

DX10F=T106*DX6F+T108*DX8F 

DX12F=D4* (DX6F+DX10F) +DD4F 

RETURN 

END 

SUBROUTINE  PER(JDR,RR) 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

COMMON/BLOCK/ AN, AM, AL,AK, PHI , T, P, KL0, IERR, X (13) 

COMMON/ TAB/DXT (12)  ,DXP(12) ,DXF(12) 

COMMON/PROP/ AVM f  R , H , U , DRT , DHT , DUT , DRP , DHP , DUP , DRF , DHF , DUF 
DIMENSION  SH (12) ,SM(12) ,HT(351) ,CT(351) 

C  DIMENSION  HT1(135) ,HT2(135) ,CTl(135) ,CT2(96) 

DATA  SfVl. 008, 16. 00, 14. 008, 2. 016, 17. 007, 28. 011, 30. 008, 31. 999, 
A  18.016,44.010,28.013,39.944/ 
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197,  4. 

,925,  5 

.668, 

6. 

427, 

I 

7. 

201, 

7.989, 

8.790, 

9.601, 

10. 

422,  11. 

,251,  12 

.087, 

12. 

930, 

J 

13. 

779, 

14.632, 

15.490, 

16.352, 

17. 

218,  18. 

,087,  18 

.958, 

19. 

833, 

K 

20. 

710, 

21.589, 

22.470, 

23.352, 

24. 

237,  25. 

,123,  26 

.011, 

26. 

901, 

L 

27. 

,791, 

28.683, 

29.577, 

30.471, 

31. 

367,  32. 

264,  33 

.161, 

34. 

060/ 

DATA  CT/5.049, 

5.029,  ! 

5.015,  5 

i.  006 

i,  4.999, 

,  4.994, 

4.990, 

A 

4. 

,987, 

4.984, 

4.982,  4 

.981,  4. 

979, 

4.979, 

9.978, 

4.978 

B 

4. 

978, 

4.979, 

4.980,  4 

.981,  4. 

984, 

4.986, 

4.990, 

4.994 

C 

4. 

999, 

5.004, 

5.010,  5 

.017,  5. 

025, 

5.033, 

5.041, 

5.050 

D 

5. 

060, 

5.070, 

5.081,  5 

.091,  4. 

968, 

4.968, 

4.968, 

4.968, 

E 

4. 

,968, 

4.968, 

4.968,  4 

.968,  4. 

968, 

4.968, 

4.968, 

4.968 

F 

4. 

,968, 

4.969, 

4.969,  4 

.969,  4. 

971, 

4.972, 

4.975, 

4.978 

G 

4. 

,982, 

4.987, 

4.993,  5 

.001,  5. 

011, 

5.022, 

5.035, 

5.050 

H 

5. 

067, 

5.086, 

5.107,  5 

.130,  5. 

156, 

5.183, 

5.213, 

7.009 

I 

7. 

,036, 

7.087, 

7.148,  7 

.219,  7. 

300, 

7.390, 

7.490, 

7.600 

J 

7. 

,720, 

7.823, 

7.921,  8 

.016,  8. 

108, 

8.195, 

8.279, 

8.358 

1 

8. 

,434, 

8.506, 

8.575,  8 

.639,  8. 

700, 

8.757, 

8.810, 

8.859 

2 

8. 

.911, 

8.962, 

9.012,  9 

.061,  9. 

110, 

9.158, 

9.205, 

9.252 

L 

9. 

,927, 

9.342, 

7.057,  7 

.090,  7. 

150, 

7.233, 

7.332, 

7.439 

M 

7. 

,549, 

7.659, 

7.766,  7 

.867,  7. 

966, 

8.053, 

8.137, 

8.214 

N 

8. 

.286, 

8.353, 

8.415,  8 

.472,  8. 

526, 

8.576, 

8.622, 

8.665 

A 

8. 

.706, 

8.744, 

8.780,  8 

.814,  8. 

846, 

8.876, 

8.905, 

8.933 

B 

8. 

.959, 

8.984, 

9.008,  9 

.031,  9. 

053, 

7.276, 

7.450, 

7.624 

C 

7. 

.786, 

7.931, 

8.057,  8 

.168,  8. 

263, 

8.346, 

8.417, 

8.480 

D 

8. 

,535, 

8.583, 

8.626,  8 

.664,  8. 

689, 

8.728, 

8.756, 

8.781 

E 

8. 

,804, 

8.825, 

8.844,  8 

.863,  8. 

879, 

8.895, 

8.910, 

8.924 

F 

8. 

.937, 

8.949, 

8.961,  8 

.973,  8. 

984, 

8.994, 

9.004, 

9.014 

G 

7. 

,466, 

7.655, 

7.832,  7 

.988,  8. 

123, 

8.238, 

8.336, 

8.419 

H 

8. 

.491, 

8.552, 

8.605,  8 

.651,  8. 

692, 

8.727, 

8.759, 

8.788 

I 

8. 

,813, 

8.837, 

8.858,  8 

.877,  8. 

895, 

8.912, 

8.927, 

0.941 

J 

8. 

.955, 

8.968, 

8.980,  8 

.991,  9. 

002, 

9.012, 

9.022, 

9.032 

K 

9, 

.041, 

9.050, 

9.058,  7 

.670,  7. 

883, 

8.063, 

8.212, 

8.336 

L 

8. 

,439, 

8.527, 

8.604,  8 

.674,  8. 

738, 

8.800, 

8.858, 

8.916 

M 

8, 

.973, 

9.029, 

9.084,  9 

.139,  9. 

194, 

9.248, 

9.301, 

9.354 

N 

9. 

.405, 

9.455, 

9.503,  9 

.551,  9. 

596, 

9.640, 

9.682, 

9.723 

0 

9, 

.762, 

9.799, 

9.835,  9 

.869,  9. 

901, 

9.932, 

8.676, 

8.954 

I'.* 
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A 

9. 

246, 

9. 

574, 

9. 

851, 

10 

.152, 

,  10.444 

t 

10.723, 

10. 

987 

,  11 

B 

11 

.462, 

,11 

.674, 

,11 

.869, 

,12 

.048, 

,12 

.214, 

12 

.366, 

12 

.505 

,12 

.634 

C 

12 

.753, 

,12 

.863, 

,12 

.965, 

,13 

.059, 

,13 

.146, 

13 

.228, 

13 

.304 

,13 

.374 

E 

13 

.441, 

,13 

.503, 

,13 

.562, 

,13 

.617, 

,13 

.669, 

13 

.718, 

13 

.764 

,13 

.808 

D 

13 

.850, 

,11 

.310, 

,11 

.846, 

,12 

.293, 

,12 

.667, 

12 

.980, 

13 

.243 

,13 

.466 

E 

13 

.656, 

,13 

.815, 

,13 

.953, 

,14 

.074, 

,14 

.177, 

14 

.269, 

14 

.352 

,14 

.424 

F 

14 

.489, 

,14 

.547, 

,14 

.600, 

,14 

.648, 

,14 

.692, 

14 

.734, 

14 

.771 

,14 

.807 

G 

14 

.841, 

,14 

.873, 

,14 

.902, 

,14 

.930, 

,14 

.956, 

14 

.982, 

15 

.006 

,15 

.030 

H 

15 

.053, 

,15 

.075, 

,15 

.097, 

,15 

.119, 

,  7 

.196, 

7 

.350, 

7 

.512 

,  7 

.670 

I 

7. 

815, 

7. 

945, 

8. 

061, 

8. 

162, 

8. 

252, 

8. 

330, 

8. 

398, 

8. 

458, 

J 

8. 

512, 

8. 

559, 

8. 

601, 

8. 

638, 

8. 

672, 

8. 

703, 

8. 

731, 

8. 

756, 

K 

8. 

779, 

8. 

800, 

8. 

820, 

8. 

838, 

8. 

855, 

8. 

871, 

8. 

886, 

8. 

900, 

L 

8. 

914, 

8. 

927, 

8. 

939, 

8. 

950, 

8. 

962, 

8. 

972, 

8. 

983, 

8. 

993/ 

C  SUBROUTINE  PER  CALLS  SUBROUTINE  EQMD  FOR  CALCULATTION  OF  MOLE 

C  FRACTIONS  AND  AND  THEIR  PARTIAL  DERATIVES 

C 


CALL  EQMD( JDR, RR) 

IF(IERR. NE.0)  RETURN 

C  SECTION  100  CALCULATES  AVM, R,H  AND  U. 

TK=T/180. 

IT=TK 

FR=TK~FLOAT ( IT) 

SH(1) =92935. 8+4. 968*T 
SH(12) =4.968*T 

C  WRITE (1, 887) SH (1) , SH(12) 

887  FORMAT ( IX , ' SH ( 1 )  SH(12) ' ,2F12.2) 

DO  105  1=2,11 
IN=35* (1-2) + (IT-5) 

SH(I)= (HT(IN) +FR* (HT(IN+1)-HT(IN) ) ) *1800. 
105  CONTINUE 

AVM=0 . 

H=0. 

DO  110  1=1,12 

IF(X (I) . LE.0. 000000001) X (I) =0.000000001 
AVM=AVM+X ( I ) *  SM ( I ) 

H=H+X(I) *SH(I) 

990  FORMAT ( IX, 'H  X  SH  AVM  HT(I)  DELT' ,6F10. 2) 

110  CONTINUE 

H=H/ AVM 

R=1.987165/AVM 


i 


V 

V 

V 


RR=R 

C  WRITE (1,*)  '**RR'  ,RR 

U=H-R*T 

IF (JDR. EQ. 0)  RETURN 
C 

DMT=DXT ( 1 ) *  SM ( 1 ) +DXT ( 1 2 ) *  SM ( 1 2) 

DHT=X (1)  *4. 968+DXT (1) *SH(1)+X(12) *4.968+DXT(12) *SH(12) 
DMP=DXP(1) *SM( 1) +DXP (12) *SM(12) 

DHP=DXP (1) *SH(1) +DXP(12) *SH(12) 

DMF=DXF ( 1 ) * SM ( 1 ) +DXF ( 1 2 ) * SM ( 1 2 ) 

DHF=DXF(1) *SH(1) +DXF(12) *SH(12) 

DO  210  1=2,11 
IN=35* ( 1-2) + ( IT-5) 

CP=CT ( IN) +FR* ( CT { IN+1) -CT ( IN) ) 

DMT=DMT+DXT ( I ) *  SM ( I ) 
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m 


DHT=DHT+X ( I ) *CP+DXT (I) *SH(I) 
DMP=DMP+DXP (I)  *SM{ I) 

DHP=DHP+DXP(I)  *SH(I) 

DMF=DMF+DXF ( I ) *SM(I) 
DHF=DHF+DXF(I)*SH(I) 

CONT INUE 
DRT— R*  DMT/  AVM 
DHT= (DHT~DMT*H) /AVM 
DUT=DHT-R-DRT*T 
DRP— R*DMP/AVM 
DHP=  (DHP-DMP*H)  /AVM 
DUP=DHP-DRP*T 
DRF— R*DMF/AVM 
DHF= (DHF~DMF*H) /AVM 
DUF=DHF-DRF*  T 

FORMAT ( IX / ' H  DHT  DELT  T',4F15.5) 

RETURN 

END 


| 


